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ABSTRACT 

The  present  work  is  an  attempt  to  put  together  the  most 
relevant  aspects  of  the  engineering  problems  involving  dis- 
tributed parameter  systems  (D.P.S.'s).   Simulation  and  optimal 
control  are  explained  in  detail  in  Chapters  II  and  III. 

The  original  contribution  of  this  thesis  is  given  in 
Chapters  V  and  VI,  where  modal  control  theory  and  a  gradient 
subroutine  that  searches  for  the  optimal  reference  coeffi- 
cients are  used.   As  a  result,  it  was  possible  to  obtain  an 
output  distribution  better  than  the  one  achievable  by  the 
known  methods.   This  technique  works  in  situations  of  strongly 
nonlinear  control  and  compensates  the  effect  of  having  the 
analyzer  and  synthesizer  approximated  by  low  order  matrices. 
It  also  makes  it  possible  to  give  higher  weight  to  some  zones 
of  the  output  distribution  in  order  to  have  a  better  local 
fit.   The  necessary  background  for  understanding  Chapters  V 
and  VI  is  given  in  Chapter  IV. 
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I.   GENERALITIES  ABOUT  THE  CONTROL  OF  DISTRIBUTED  PARAMETER 
SYSTEMS 


A.   HISTORY 

Distributed  parameter  systems  have  existed  in  Nature 
since  the  beginning,  and  many  of  these  systems  were  natu- 
rally stable. 

For  example,  during  the  period  of  pre-living  beings  dra- 
matic changes  in  temperature  originated  no  less  dramatic 
variations  of  the  lithosphere.   These  earth's  convolutions, 
together  with  many  subsidiary  effects,  tried  then  to  make 
changes  in  different  places  until  a  practical  equilibrium  was 
obtained.   What  was  this,  other  than  a  distributed  parameter 
system? 

With  the  advent  of  life  in  this  planet  other  D.P.S.'s  ap- 
peared, such  as  the  nitrogen  cycle,  in  which  different  forms 
of  life  receive  nitrogen  and  give  it  up  to  other  forms  which, 
in  turn,  originate  new  nitrogen. 

In  the  beginning  of  this  century  much  attention  was  dedi- 
cated to  the  mathematical  description  of  D.P.S.'s  such  as 
transmission  lines.   With  the  entrance  on  the  scene  of  modern 
control  theory  and  integrated  circuits  concepts,  some  attempts 
were  made  toward  a  generalization  of  those  ideas,  but  practi- 
cally until  the  1950' s  not  much  had  been  realized.   By  that 
time  a  certain  number  of  systems  with  time  delays  could  be 
analyzed  using  what  are  called  conventional  control  techniques, 
but  the  scope  of  these  techniques  was  quite  limited. 
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In  the  early  1960 's  some  technical  papers  were  published 
in  an  attempt  to  apply  the  concepts  of  optimal  control  to 
D.P.S.'s  and,  from  then  on,  a  highly  geometric  rate  of  tech- 
nical paper  production  took  place. 

A  new  kind  of  approach  to  the  control  of  D.P.S.'s  was 
carried  on  by  Murray-Lasso  [Ref.  45]  and  Gould  [Ref.  26]  in 
1965,  using  the  generalization  of  the  concepts  of  modal  con- 
trol introduced  by  Rosenbrock  [Ref.  53] .   This  approach  is 
treated  in  depth  in  Chapters  IV,  V  and  VI.   References  62  and 
18  also  contain  very  useful  information  in  this  field. 

In  196  8  the  Int.  Journal  of  Control  published  an  exten- 
sive bibliography  by  Wang  [Ref.  72],  one  of  the  scientists  in 
the  U.S.A.  who  has  contributed  most  significantly  to  the  de- 
velopment of  the  D.P.S.'s  theory.   In  November  1969,  the 
Aerospace  Research  Laboratories  published  the  "Survey  of  op- 
timal Control  of  Distributed  Parameter  Systems,"  by  A.  C. 
Robinson  [Ref.  52],  which  contains  261  references  and  an  anal- 
ysis of  various  ways  of  implementing  the  optimal  control  of 
such  systems.   Special  reference  ought  to  be  made  to  the 
contributions  of  Lions  [Ref.  40]  and  Butkovskiy  [Ref.  5], 
which  respectively  in  France  and  U.S.S.R.  established  founda- 
tions of  a  scientific  treatment  of  the  optimal  control  of 
D.P.S. 's. 

In  1969  and  1970  numerous  technical  papers  on  D.P.S.'s 
were  published.   Many  of  them  may  be  found  in  the  "Proceedings 
of  the  Joint  Automatic  Control  Conference,"  "Int.  Journal  of 
Control"  and  "Simulation." 
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Almost  all  the  above  references,  although  necessary,  are 
very  theoretical  and  not  much  has  been  done  by  the  control 
engineer  in  order  to  implement  those  methods  or  to  describe 
them  in  an  easy  language.   It  is  a  purpose  of  this  thesis  to 
fill  a  little  of  this  existing  gap. 

History  is  not  only  the  study  of  the  past  events  but  also 
the  careful  analysis  of  them  in  order  to  extrapolate  for  the 
future.   It  seems  reasonable,  on  the  basis  of  past  history  to 
make  some  predictions  as  to  future  studies  of  D.P.S.'s. 

i)   From  a  technological  point  of  view  it  is  a  known  fact 
that  the  impact  of  control  science  in  the  improvements 
accomplished  in  the  last  two  decades  constituted  a  "sine 
qua  non"  factor  for  these  trends.   At  the  actual  rates 
of  development  of  new  D.P.S.'s  (as  mentioned  before  the 
integrated  circuits  are  a  characteristic  example)  and 
of  the  capability  to  describe  and  control  them,  it  seems 
that  within  a  few  years  the  control  science  will  be  ap- 
plicable to  almost  every  different  aspect  of  the  existing 
technology . 
ii)   The  actual  and  very  important  ecological  problems,  namely 
water  and  air  pollution,  if  mathematically  described  will 
assume  the  form  of  partial  differential  equations.   The 
consequences  of  being  able  to  deal  with  such  problems  in 
a  scientific  way  can  be  easily  implied. 
iii)   It  seems  realizable  that  partial  differential  equations 
may  be  derived  that  approximate  the  behavior  of  some 
segments  of  society,  at  least  during  certain  periods  of 
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time.   When  such  equations  'are  found  it  should  be  pos- 
sible to  use  computer  optimization  (parameter  identifi- 
cation) to  determine  values  or  expressions  for  the 
coefficients.   These  coefficients  are  but  the  histori- 
cal constants  the  historians  have  been  looking  for. 
iv)   Man  is  a  complex  being,  and  his  free  nature  leads  to 

impulsive  behavior  at  random  intervals.   The  effect  of 
such  behavior  on  society  is  often  observable.   It  is 
reasonable  to  think  that  the  phenomena  of  man's  behavior 
may  in  some  way  follow  a  Gaussian  distribution,  as  do 
many  natural  phenomena.   If  equations  can  be  obtained  to 
represent  society  as  a  D.P.S.,  then  it  may  be  possible 
to  treat  the  impulsive  aspects  of  man's  behavior  as 
Gaussian  noise  and  perhaps  the  consequences  of  such  be- 
havior can  be  minimized  by  building  into  the  social  sys- 
tem something  similar  to  a  Kalman  filter. 
v)   As  progress  is  made  in  developing  equations  for  various 
sub-systems  within  the  social  systems,  it  may  well  be 
that  some  sub-systems  may  be  found  uncontrollable,  others 
controllable.   Then  the  theory  of  D.P.S.'s  may  well  be- 
come an  important  tool  for  redesign  of  social  systems. 
Finally,  in  order  to  give  a  geometrical  interpretation 
of  the  behavior  of  a  D.P.S.,  Figs.  1.1  and  1.2  are  included, 
which  show  how  each  single  input  contributes  to  the  value  of 
the  output  at  every  point. 
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Figure  1.1.   A  D.P.S.:   Heating  of  a  Rod 
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Figure  1.2.   A  More  General  Model 

Note:   The  output  transformation  uses  the  values  measured  by 

some  existent  sensors  and  converts  them  into  the  desired 
output. 
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B.   HOW  TO  CONTROL  DISTRIBUTED  PARAMETER  SYSTEMS 

The  ordinary  differential  equations  used  by  the  control 
engineer  to  describe  lumped  parameter  systems,  with  time  as 
independent  variable,  did  not  prove  sufficiently  accurate  to 
deal  with  complex  systems  where  either  the  input,  or  the  out- 
put, or  both,  are  also  function  of  other  parameters,  generally 
of  spatial  nature.   Therefore,  a  new  theory  was  developed  in 
order  to  extend  some  of  the  ideas  so  well  understood  when 
dealing  with  ordinary  differential  equations  to  the  much  more 
difficult  problems  described  by  partial  differential  equations 
and  integral  equations. 

There  are  four  basic  approaches  to  the  control  problem  of 
D-.P.S.'s;  however,  none  of  them  is  a  general  method  applicable 
to  all  types  of  these  systems: 

1.  Use    of   Time-Delay   Techniques 

The    use    of   time-delay   techniques    [Ref.    64]    has   been 
common   when   dealing  with   problems    such   as   those    that    arise    in 
long  pneumatic   control    lines,    heat   transfer   in   nuclear   plants, 
production   delay    in    assembly    lines    and   in   general    to  every 
kind   of    control    systems   where    time-delays    are    involved.      The 
conventional    analysis    and   design    techniques    applied   to   these 
systems    requires    too  much   work    and    for   this    reason    they    are 
now   studied    almost   entirely   by    computer   simulation. 

2 .  Reduction    of   the    Partial    Differential    Equations    to 
Ordinary    Differential    or   Difference    Equations 

Once    the    ordinary   differential    or   difference   equations 

are    obtained    lumped  parameter   techniques    must   be    used.       From   a 
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control  point  of  view  the  most  representative- methods  used  to 
obtain  the  above  equations  are  as  follows: 

a.  Space  Quantization  (DSCT) 

The  space  derivatives  are  replaced  by  finite  dif- 
ferences but  the  time  derivatives  are  maintained.   One  of  the 
biggest  advantages  of  this  method  is  that  it  permits  an  analog 
modeling  of  the  system.   The  accuracy  can  be  improved  by  in- 
creasing the  number  of  sections,  which  must  be  adequately 
isolated.   When  the  number  of  sections  becomes  too  large  for 
the  required  accuracy,  another  method  must  be  chosen. 

b.  Time  and  Space  Quantization  (DSDT) 

Use  is  made  of  the  well-known  techniques  for  the 
numerical  solution  of  partial  differential  equations.   Al- 
though quite  valuable  from  an  analytic  point  of  view  and  for 
simulation,  it  seems  a  little  lengthy  if  one  is  trying  to  ap- 
ply it  to  control  design. 

c.  Laplace  Transform 

The  Laplace  transform  method  uses  the  fact  that 

Tr3f(x,t),    3F(x,s)    ,  Tr8F(x,t)1     _ ,    .    cl       n+x    . 
L[- —  '   ]  =  \'         and  L[ ~         ]    =  sF(x,s)  -  f (x,0  )  ,  when- 

d  X  a  X  du 

ever  f(x,t)  is  transformable,  L[f(x,t)]  =  F(x,s)  and  — t^t 

exists. 

The  problems  that  can  be  treated  by  Laplace  trans- 
formations are  quite  numerous,  the  principal  objection  being 
the  intricate  expressions  to  which  one  arrives  when  dealing 
with  systems  somehow  more  complex  (i.e.  multidimensional 
systems) . 
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3.  Modal  Control 

The  modal  control  approach,  Refs.  26  and  45,  consists 
in  replacing  the  partial  differential  operator  with  an  in- 
finite order  matrix  which  is  approximated  by  another  matrix 
of  order  N<°°.   Through  adequate  matrix  manipulations  it  is 
possible  to  uncouple  the  interconnected  inputs  and  outputs, 
such  that  separate  loops  are  established,  each  of  which  can  be 
treated  by  lumped  parameter  control  techniques. 

Its  most  noteworthy  advantage  is  the  great  insight  that 
it  can  give  to  a  problem.   The  range  of  operation  is,  however, 

limited  because  it  can  be  applied  only  to  completely  continuous, 

2 

bounded  and  linear  or  linearized  operators.    Also,  sometimes, 

the  truncation  is  problematic,  either  because  some  of  the 
higher  order  eigenvalues  may  be  of  fundamental  importance  or 
because  the  convergence  of  these  eigenvalues  may  be  too  slow. 
In  this  case  the  system  would  require  a  high  order  matrix 
representation,  which  is  not  very  acceptable  in  practice. 
The  last  drawback  seems  to  be  the  difficulty  that  often  arises 
in  the  computation  of  the  eigenvalues  and  eigen functions.   A 
detailed  analysis  of  this  theory  is  given  in  Chapter  IV. 

4 .  Optimal  Control 

Optimal  control  is,  by  far,  the  technique  about  which 
the  largest  amount  of  work  has  been  published.  This  comes  from 
the  fact  that  the  control  of  systems  described  by  partial  dif- 
ferential equations  varies  remarkably  with  changes  in  the 
initial  and  boundary  conditions,  with  the  definition  of  cost 


2  ... 

See  Appendix  B  for  definitions. 
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function  and  from  system  to  system.   Some  more  or  less  gener- 
alized solutions  have  been  proposed,  normally  according  to  the 
classification  of  the  system  operator  within  the  three  fol- 
lowing categories:   elliptic,  parabolic  or  hyperbolic  partial 
differential  operator.   Optimal  control  will  be  discussed  in 
a  certain  depth  in  Chapter  III;  it  seems  to  be  the  approach 
of  greatest  future.   This  is  due  to  the  increasing  knowledge 
on  how  to  describe  the  cost  function,  giving  convenient 
weight  to  each  of  its  components,  and  also  due  to  the  tendency 
for  the  control  systems  to  become  highly  computerized  because 
of  the  low  cost  of  integrated  circuitry. 

The  four  basic  approaches  just  described  are  not  com- 
pletely independent  of  each  other.   Effectively,  in  many 
circumstances,  all  of  them  are  ultimately  reduced  to  the  solu- 
tion of  problems  involving  ordinary  differential  or  difference 
equations.   There  is,  however,  an  important  reason  for  the 
chosen  criterion,  which  can  be  better  explained  giving  the 
following  example.   Consider  para.  2;  it  is  obvious  that  after 
the  transformation  is  accomplished  the  lumped  parameter  optimal 
control  techniques  can  be  used;  the  difference  from  para.  5  is 
in  the  way  the  optimal  control  is  derived.   In  the  first  case 
the  optimal  control  law  is  derived  from  an  ordinary  differen- 
tial or  difference  equation,  while  in  the  last  one  it  is 
derived  directly  from  the  given  partial  differential  equation. 
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II.   SIMULATION  METHODS  FOR  DISTRIBUTED  PARAMETER  SYSTEMS 

A.   INTRODUCTION 

The  distributed  parameter  part  of  a  physical  system  is 
most  accurately  described  by  partial  differential  equations. 
Analysis  and  design  of  such  systems  requires  solution  of  the 
partial  differential  equations.   Solution  by  analytical  meth- 
ods can  be  quite  laborious,  and  such  solutions  are  not  help- 
ful in  design  problems;  indeed  they  may  not  be  convenient  for 
many  system  analysis  problems. 

An  alternate  approach  is  the  simulation  of  the  distributed 
parameter  system  in  a  computer  (analog,  digital  or  hybrid)  and 
the  use  of  this  simulation  model  in  computer  studies  planned 
for  analysis  and/or  design.   A  number  of  different  simulation 
techniques  have  been  developed  for  this  purpose.   Some  of  the 
methods  have  been  chosen  to  fit  a  specific  type  of  computing 
facility,  others  have  been  designed  to  solve  a  specific  class 
of  partial  differential  equations,  and  still  others  are  general 
methods  that  can  be  applied  to  a  variety  of  problems. 

The  purpose  of  this  chapter  is  to  classify  and  list  most 
of  the  available  simulation  methods,  giving  also  some  of  the 
important  facts  associated  with  each  of  them.   No  detailed 
comparison  or  evaluation  is  attempted  since  the  choice  of  a 
method  is  guided  largely  by  available  computer  facilities, 
the  nature  of  the  specific  problem  to  be  solved,  the  accuracy 
needed  in  the  solution,  and  the  constraints  on  cost  and  time. 
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However  the  chapter  is  concluded  with  a  summarizing  table 
which  indicates  some  of  the  pertinent  considerations  associ- 
ated with  each  method.   This  table  could  be  used  as  a  start 
of  a  critical  comparison  and  evaluation. 

B.   A  GENERALIZED  CLASSIFICATION 

It  can  be  said  that,  in  general,  the  partial  differential 
equation  describing  a  distributed  parameter  system  is  part  of 
an  overall  control  system.   In  order  to  analyze  the  behavior 
of  this  complete  system  it  is  necessary  to  know  how  to  solve 
the  corresponding  p.d.e.   This  can  be  done  either  by  analyt- 
ical methods  or  by  simulation. 

Simulation  has  the  great  advantage  that  once  the  system  is 
represented  for  a  specific  input  and  initial  conditions  any 
change  of  conditions  is  easily  realizable;  such  is  not  the 
case  for  the  analytical  computation. 

The  techniques  for  dealing  with  multidimensional  systems 
are  only  extensions  of  those  used  for  unidimensional  ones. 
Due  to  the  fact  that  the  simulation  of  those  systems  involves 
an  enormous  amount  of  hardware,  if  implemented  by  analog  meth- 
ods, most  multidimensional  simulations  have  been  done  in  the 
digital  computer  and  sometimes  in  the  hybrid  computer.   They 
are,  however,  tremendously  time  consuming  and  for  this  reason 
only  the  two-dimensional  systems  can  be  modeled  with  reason- 
able, if  not  great  accuracy  [Ref.  49].   The  objective  of  a 
simulation  will  be  to  have  for  specific  positions  the  output 
as  a  function  of  time  or,  for  each  instant  of  time,  the  output 
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as  a  function  of  the  coordinates.   In  a  two-dimensional  sys- 
tem the  solution  is  normally  obtained  keeping  two  coordinates 
fixed,  changing  the  other  one  and  repeating  this  for  several 
combinations  of  points  and  time  instants  until  a  complete  grid 
is  obtained.   Actually  some  oscilloscopes  are  already  capable 
of  providing  multidimensional  pictures.   This  thesis  will  deal 
only  with  one -dimensional  systems,  for  which  the  following 
methods  of  solution  may  be  considered: 

DSDT  -  discrete-space,  discrete-time 

DSCT  -  discrete-space,  continuous -time 

CSDT  -  continuous-space,  discrete-time 

TSCT  -  transformed-space,  continuous -time 

TSDT  -  transformed-space,  discrete-time 

Schuchmann  [Ref.  59]  gives  a  brief  and  clear  analysis  of 
these  methods,  involving  error  and  stability  considerations, 
as  well  as  the  relative  advantages  of  each  one.   He  doesn't 
consider  the  three  possibilities  resulting  from  taking  the 
time  Laplace  transform,  explaining  that  this  would  imply  the 
need  for  getting  its  inverse;  this,  although  feasible  is  very 
time  consuming.   However,  when  considering  the  solution  at  few 
points  (sometimes  only  one)  the  technique  seems  to  be  good,  in 
particular  if  using  infinite  product  expansions  [Ref.  22]. 
Such  technique  gives  the  solution  as  a  product  of  terms 
(truncated  infinite  series)  in  "s",  each  of  them  easily  im- 
plementable  in  an  analog  computer;  it  will  be  discussed  in  more 
detail  in  Section  F  of  this  chapter. 

1.   DSDT  -  Discrete-Space,  Discrete-Time 

With  the  very  fast  digital  computers  actually  exis- 
tent and  also  because  of  their  enormous  accuracy,  popularity 
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and  availability,  most  of  the  work  has  the  tendency  to  be 
done  by  the  DSDT  method.   Also  with  the  actual  video-display 
units  it  is  possible  to  get  CRT  pictures  and  change  easily 
the  parameters  without  need  for  an  analog  or  hybrid  computer. 

The  basic  concepts  of  the  discretization  are  very 
simple  and,  because  of  its  importance,  all  of  the  next  sec- 
tion will  be  devoted  to  this  problem. 

2.   DSCT  -  Discrete-Space,  Continuous -Time 

This  is  the  method  generally  used  when  there  is  avail- 
able only  the  analog  computer  and  the  transformation  tech- 
niques described  previously  in  this  chapter  are  difficult  to 
implement.   If  the  number  of  nonlinearities  or  the  required 

accuracy  is  large,  the  lack  of  enough  multipliers  or  other 

3 

components  m  the  computer  does  not  permit  solution.    To 

avoid  such  situation  the  multiplex  method  has  been  tried;  it 
consists  in  switching  on  and  off  sequentially,  with  digital 
signals,  the  parallel  operation  of  the  circuits  corresponding 
to  each  of  the  discretized  positions  and  in  using  the  same 
analog  hardware  to  simulate  the  nonlinearities  of  each  branch 
of  the  parallel  combination. 

In  Fig.  2.1  is  shown  how  the  heat  equation 

oX 


Some  companies  work  with  huge  analog  computers,  as  for 
example  English  Electronics  Co.  which  operates  the  1,500- 
amplifier  Saturn  computer. 
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can  be  simulated  in  the  analog  computer.   C  is  the  heat  ca- 
pacity and  the  heat  flow  F  is  given  by 

F=-k|H  (2.2) 

If  the  temperature  "u"  is  measured  at  integer  stations  and 
the  flux  at  half -integer  stations,  Eq.  2.2  may  be  written  in 

u  -u  _, 

discrete  form  as  F   ,  ,_  =  -  k   ,  ,„   n.  n and  Eq.  2.1  will 

n-i/2       n-1/2    Ax  ^ 

^^   9  ,v9u.|      9F,       Fn+l/2~Fn-l/2 

become  ~ — (k-^ — )   =  — ^ —   =  -  —. — 

dx      dx  |n     9x'n  Ax 

from  which 

r      ^J}    Fn+l/2"Fn-l/2  .„  _. 

n  dt  =        Ax  U-J' 

The  above  procedure  may  be  easily  implemented  accord- 
ing to  the  schematic  diagram  on  the  following  page. 
3.   CSDT  -  Continuous-Space ,  Discrete-Time 

This  technique  is  advisable  for  the  simulation  of  the 
flow  equation 

8t  +  V3x"=  f(x't} 
which,  in  discrete  form,  is  written  as 

du  (x)       ,  f  (x) 

— -5 = rrlu  (x)-u   ,  (x)  ]  +  — (2.4) 

dx       vAt   n     n-1        v 

From  this  equation  it  is  obvious  that  At  can  be  de- 
creased without  increasing  the  number  of  components,  the  only 
limitation  being  the  stability  requirement  (l/vAt<l) .   The 
analog  computer  with  adequate  number  of  track -hold  units  or 
the  hybrid  computer  are  the  most  efficient  ways  of  solution  of 

this  problem. 
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Figure  2.1 


Analog  Computer  Simulation  of  the  One 
Dimensional  Heat  Equation 
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4.  TSCT    -   Transformed-Space,    Continuous -Time 

TSCT   is    a  method   applicable    to   systems    represented  by 
self-adjoint   operators    (i.e.    heat    and  wave    equations)    and, 
through    convenient   assumptions    also   to  non-self    adjoint    sys- 
tems.     The    concepts    involved   are    intimately    related  with   those 
of   the   modal    control   which  will   be    studied   in    depth   in    Chap- 
ters   IV   and   V.      The    most    remarkable   example    of   this    method  is 
the   Bubnov-Galerkin   transformation   described   in   Appendix  A. 
Another   representative  example    is    the    one    developed   and  ex- 
plained  in   very    comprehensive    terms   by  Hong   and    Larson    [Ref. 
31]    which    can   be    applied   to    any   set   of    first-order  partial 
differential   equations    in   normal    form.       It    gives    the    solution 
as    a   continuous    function   of   distance    for   any   value    of   time. 
The    only   trade-off   is    the    requirement    for   the    computation    of 

the    inverse    transform.      However,    because    the    output    is 

2 
U (t, s) =f n~f , s+1/2 I (f„s    )-    ...,    where    f    (n=0 ,1 ,2 , . . . )    are   ex- 
clusively   functions    of   time,    the    inverse    transform   is    very 
easily    obtained.       In   the   example    given    in   the   paper   the    analog 
solution   of    a   single   nonlinear    first-order   p.d.e.    required 
only   2    integrators,    3   multipliers,    5    amplifiers    and   1    divider. 
Because    of   the  excellent    features    of   this   method    a  more    de- 
tailed  analysis   will   be    given    in    Section   E. 

5 .  TSDT    -   Transformed-Space    Discrete-Time 

TSDT  is  an  extension  of  the  above  method  and  it  is 
becoming  increasingly  more  popular.   One  example  of  such  a 
method  is  implemented  in  the  computer  run  No.  7  of  this 
thesis. 


30 


C.   DIGITAL  COMPUTER  SOLUTION  OF  PARTIAL  DIFFERENTIAL 
EQUATIONS  (DSDT) 

The  principal  source  of  information  for  this  section  is 
the  work,  "Numerical  Solution  of  Partial  Differential  Equa- 
tions/' by  G.  D.  Smith  [Ref.  63],  an  excellent  book  written 
in  very  comprehensive  terms.   For  sophisticated  problems  more 
advanced  texts  are  recommended,  such  as  Refs.  1  and  17.   Due 
to  their  relevancy  only  second-order  p.d.e.'s  will  be 

considered. 

4 
1.   Parabolic  Equation 

2 

Consider    the   parabolic  p.d.e.    ■%=■  =   K — =-  which 

2  X  U 

normalized   using    the    transformation   t=kT/L    ,    x=—  and  u=tt— , 

L       UQ 

where  Un  is  some  given  value  of  U  (generally  maximum  or  mini- 
mum) at  zero  time.   The  resulting  equation  is 


8U  _  d2u 

at  '  a  2 

X 


(2.5) 


Written    in    discrete    form   the    above    equation    becomes 

Ui,j+1        Ui+l,j    -   2ui,1    +    ui-l,j 
k  ~  h2 


(2.6) 


Given    the    second   order  quasi-linear  p.d.e. 

a~ — =-  +   Dt: — k —  +    c — k-  +  e=0   where    a,b,c,e    may   not    be    functions    of 

x  °y 

second   or   higher   order   derivatives    (according   to   the    defini- 
tion   of   quasi-linearity) ,    if   b2-4ac    is    equal    to    zero   the 
equation    is    said   to   be   parabolic,    if    greater   than    zero   is 
hyperbolic    and    is    elliptic  when    less    than    zero. 
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or 


u.  ..-,  «  u.  .  +  r(u.  ,  .  -  2u.  ,  +  u. .,  •) 


(2.7) 


Geometrically  the  interpretation  of  this  formula  is  simple, 
as  can  be  seen  from  Fig.  2.2. 
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Figure  2.2    Grid  Solution  of  a  Parabolic  Differential 
Equation 

Here,  all  the  values  at  t=0  and  x=0  or  I    are  known.   The 
computation  is,  therefore,  straightforward  and  does  not  re- 
quire iterations.   In  the  case  of  mixed  boundary  conditions 

of  the  type  au+b-^=c,  where  n  is  the  outward  normal  to  the 
ir       8n 

surface  and  a,b,c  are  functions  of  the  space  coordinates  and 
time,  only  the  values  of  u  at  t=0  are  known.   In  this  case, 
the  solution  of  the  problem  becomes  only  slightly  more  dif- 
ficult, since  it  is  necessary  to  evaluate  with  a  simple 
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formula  the  boundary  at  x=0  for  all  discrete  times,  at  the 
start  of  the  computation  of  the  elements  of  each  row  [Ref. 
47]  . 

This  may  be  done  introducing  a  new  difference  equation 
based  on  the  boundary  condition,  as  for  example 

Ul'j6x  Vj  =  h(u0,j-V)  °r  "1,J26xU"X'J  =  h(u0,j-V»  •  which 
is  more  accurate. 

The  last  equation  has  u  ,    as  unknown  and  therefore 

there  is  need  for  still  another  one  assuming  that  the  heat 

equation  is  satisfied  at  x=0. 

The  explicit  method,  although  the  simplest  one,  has  the 

2 

inconvenient  of  only  being  convergent  for  0<k/h  £1/2.   For 

this  reason  some  other  methods  were  developed  which  are  un- 

2 
conditionally  convergent  for  all  values  of  r=k/h  . 

The  most  popular  one  is  the  Crank -Nicolson  implicit  meth- 
od, which  is  represented  by  the  equation: 

Ui,i  +  l-Ui,i.lfUi+l,j  +  l-2Ui,j  +  l  +Ui-l,jH-l 
k        2*  h2 


u. , ,  . -2u.  .+u.  .  . 

2 
h 


(2.8) 


or 


-r 


u.  ,   ^n  +  (2+2r)u.  .,,-ru.,,  .,n=ru.  ,  .+  (2-2r)u.  .  +  ru.  .  . 
1-1,3+1        i/3+l    1+1,3+1    i-l/D         !/3    i+l/D 

(2.9) 
and  assuming  that  the  system  is  discretized  in  N-l  sections, 
for  each  instant  of  time,  N-l  simultaneous  equations  must  be 
solved.  The  best  way  of  solving  this  system  of  equations  is 
using  any  of  the  well-known  iterative  techniques  applied  to 
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tne    digital    computer.      The   best  known   techniques    are    Jacobi, 
Gauss-Seidel    and   successive    over-relaxation    (s.o.r.).      The 
first  one    is    very   slow   and   for   this    reason,    in   general,    only 
the    other   two   are    used,    rewriting  Eq.    2.9    as 


u.     .  . ,    =  u.     .    +   ttt{  (u.     ,       Ll  -2u.     . , ,    +   u.  ,  .     .,,) 
1,3+1           l,^         2           1-1,3+1  i,D+l  1+3,3+1' 

+    (u.     .     .  -    2u.      .    +    u.^,     .)  } 
i-l,3  i/D  1+1.3 


(2.10) 


and  dropping  the  index  j+1,  Eq.  2.10  may  now  be  expressed  as 


u.  =  =-r(u.  ,  -  2u.  +  u.  ,,)  +b. 
i    2    i-l      i     i+I     l 


whe  re  b.  =  u.  .  +  Tr(u.  ,    -2u.    +  u .  ,  ,  . ) 
i     1,3    2    1-1/3      if]     i+l,: 


(2.11) 
(2.12) 


According  to  the  above  notations  and  also  denoting  the  itera- 
tions by  superscripts,  the  Gauss-Seidel  scheme  is  obtained: 

b. 


u 


(n+1)  _    r 


2(l+r) l  i-l 


{u.^t1*  +  ufn),}  + 


i+1    1+r 


It   converges    for    all    r>0.       In   matrix    form,    making    p   = 


(2.13) 


2  (1+r) 


(n+1) 
1 

0 

p 

o 

'  4n)  1 

2 

= 

0 

2 

P 

P 

(n) 
u* 
2 

(n+1) 
Q3 

0 

3 

P 

p2       P  o 

(n) 
U3 

'(n+1) 
Tl-1 

0 

N-l 
P 

N-2 

P 

1  (n) 

+  c 


(2.14) 


The    successive    over- re  lax at ion    scheme    is    the    fastest    one    and 
is    obtained   directly    from  Eq.    2.13   by    adding    and   subtracting 

u. (n)     to   the    right-hand    side. 

The    following  expression    is    obtained: 


u.(n+1)=u.(n)+U[ 
1  1 


b. 


-    u(n)] 
2(l+r)L"i-l       ' "i+lJ"(l+r)  i      J 


-{u^+uf^H^ 


(2.15) 


where    oo   is    the    relaxation    factor    (usually    l<to<2)  . 
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For  maximum  rate  of  convergence 

%  =   "ZZZZZZZ     where   H   =   xfr   cos   ff  (2-!6) 

i+/d-y2) 

2 .   Hyperbolic  Equations 

The  most  commonly  used  method  for  the  numerical  solu- 
tion of  hyperbolic  equations  is  the  method  of  "characteristics" 
and  it  consists  in  integrating  the  given  p.d.e.  along  chosen 
directions  (when  dealing  with  two-dimensional  systems;  if 
three  dimensional  systems  are  considered  the  integration  is 
in  chosen  planes)  for  which  the  partial  derivatives  are  re- 
duced to  ordinary  derivatives. 

The  mathematical  explanation  is  as  follows:   given 
the  quasi-linear  equation 

a§+brif7+cl7  +  e  =  0  (2-17) 

2  2 

with  b    -4ac>0    and   defining   —  =  p,r—  =   q,       g— j  =    r»    JTxTy"  = 

.2  X 

8    u 
and   - — =■  =   t    ,     can   now   be   written: 

dp   =   ZE.  dx   +   |£  dy   =   rdx   +    sdy 
^         3x  8Y 

and 

dq   =    a£L  dx   +    23.  dy   =    sdx   +    tdy  (2.18) 

8*  8Y 

Also    from  Eq.    2.17:       ar+bs+ct+e   =    0,    which   is    used   to- 
gether  with   Eq.    2.18    to  eliminate    r    and   t,    giving    the    final 
expression 

s{a(dY)2    _   b(dY)    +    c}    -    {A  dY   +      dq   +   edy}    =    0  (2.19) 

1       dx  dx  dx   dx  dx  dx 
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dv 
Choosing  -g-   such  that  the  polynomial  inside  the  left  {  } 


dx 
is  zero,  it  will  be  necessary  to  solve  only 

dP  dy    dq    dy  = 

dx  dx     dx     dx  v     ' 

where  -v^-  is  known  and  takes  two  real  values. 

A  detailed  procedure  for  the  computer  implementation  of 
this  technique  is  given  in  Ref.  63,  pp.  101-102. 

The  method  of  characteristics  is  the  most  accurate  one 
for  solving  hyperbolic  p.d.e.'s  and  whenever  discontinuities 
are  involved  is  the  one  that  must  be  used.   However,  if  no 
discontinuities  are  present,  a  convergent  finite , differences 
method  is  much  easier  to  implement  and  will  give  generally 

adequate  results.   One  scheme  that  may  be  used  when  dealing 

.2     .2 

with  the  equation  » — T  =  — s-  is 

a     Z  r\^Z 

X      8t 

u.       L-,=r2    u.     ,         +2(l-r2)u.     .    +r2    u.^,     .    -u.     .     ,  (2.21) 

i,3+l  i-l,D  i,D  1+1/3         i/D-1 

6x 
with  r=-p-r-,  which  is  convergent  for  r<l. 

3.   Elliptic  Equations 

Laplace's  and  Poisson's  equations  are  the  two  best 
known  equations  of  the  elliptical  type.   These  equations  have 
the  peculiarity  of  being  always  integrable  in  a  closed  area 
(if  a  two-dimensional  equation)  or  volume  (if  a  three- 
dimensional  equation)  on  which  boundaries  the  function  or  its 
gradient  are  known. 

For  example  Poisson's  equation 

9U  +  dJJ   =  f{XfY)  (2#22) 

ax    8YZ 
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X  Y 

may   be    represented   in   normalized    form  making   x=—   ,    y=—   , 

Li  L 

and   u=— j   .       One   possible    discrete -vers ion   may  be 
.  IT 


u-:_i.i    a+ua     i     ^+u-     •J.i+U-      •    i-4u.     •    =   h    f(x.y.) 


(2.23) 


which  must  be  verified  for  every  point  inside  the  boundary. 
In  general,  if  there  are  N  points  inside  the  boundary  a  set  of 
N  equations  will  result,  linear  or  nonlinear,  in  the  same  way 
as  the  original  p.d.e. 

When  the  boundaries  are  curvilinear  and  the  function  val- 
ues or  its  derivatives  cannot  be  represented  accurately  by 
the  chosen  mesh,  either  a  finer  mesh  must  be  used  or  the  fol- 
lowing equation: 


2u. 


2u 


+ 


B 


2u 


+ 


07(1+07)         02(l+02)         1+01 


+ 


2u4  11  2 

IT07-   2<07+   ^)uo  =  h    fo 


(2.24) 


according  to  Fig.  2.3,  where  f   is  the  value  of  f  at  (0,0). 


(0<@,<1 ) 
(0<©,<1) 


Figure  2.3.  Grid  for  an  Elliptic  p.d.e.  with  Curvilinear 
Boundaries 
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D.   APPLICATION  OF  THE  CONCEPT  OF  "CHARACTERISTICS"  IN  HYBRID 
COMPUTERS  (DSCT  OR  CSDT) 

There  exists  many  different  ways  of  solving  p.d.e.'s  in 
a  hybrid  computer.   Maybe  the  most  representative  of  these 
procedures  is  the  extension  of  the  concept  of  the  character- 
istics described  in  the  last  section.   Such  an  extension  was 
originated  by  Vichnevetsky  and  Associates  [Refs.  69  and  70]. 
It  is  explained  as  follows:   given  the  p.d.e. 

|H  +  f(X/t)|^=  g(u,x,t)  (2.25) 

dx 
with  given  initial  and  boundary  conditions  and  making  -rr-  = 

f(x,t)  it  is  possible  to  write 

du    3u  ,  3u    dx    3u  .  ^du  ,~  orv 

3-  =  ttt-  +  77—  *  jt  =  -kit   +  f  t~  (2.26) 

dt    9t    8x    dt    9t     dx 

Therefore 

^  =  g(u,x,t)  (2.27) 

which  is  an  ordinary  differential  equation,  integrable  by 
analog  computer. 

For  the  boundary  condition  u(0,t),  Eq.  2.27  can  be  inte- 
grated at  high  speed  while  the  integration  of  x  is  done 
slowly  (characteristic  line) .   Geometrically,  for  L  the 
length  of  the  system  this  method  is  exemplified  in  Fig.  2.4. 
If  the  highest  integration  speed  is  much  greater  than  the 
lower  one  the  time  can  be  taken  as  a  constant  for  each  fast 
sweep;  otherwise  a  correction  can  be  implemented.   For  fur- 
ther information  a  specific  case  is  worked  out  with  great 
detail  in  Ref.  70.  - 
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X  , 


Figure    2.4.       Geometrical    Representation    of   the    Method 
of    Characteristics    in   Hybrid   Computer 


E.       ANALOG    SIMULATION    OF    UNI DIMENSIONAL    D.P.S.'s    BY 

LAPLACE    TRANSFORMING    THE    SPACE    COORDINATES     (TSCT) 

Following    the    same    procedure    as    used  by   Hong   and   Larson 

in    Ref.    31,    consider   the    system   defined   by 


|H+    f(t)|H+g(t)u-0 


(2.28) 

with  initial  and  boundary  conditions 

u(0,x)  =  ae"g(o)x  (2.29) 

u(t,0)  =  a 

Normalizing  the  length  and  by  defining  u  (t,l)=b  the  distance 

Laplace  transform  of  u  will  be 


-sx 


U  (t,s)  =  /  e    udx  =  / 


-sx 


udx 


(2.30) 
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and,  after  substituting  in  Eq.  2.28  it  comes  out 


dD'£'S)  +  f  (t)[/  e"SX  dU^'s)dx)+g(t)  DCs)  =  0       (2.31) 

o 

or,  because  the  integral  is  identical  to  zero  for  x  greater 

than  1 : 

dU^'S)  +  (g  (t)  +  sf  (t))  U(t,s)  =  f(t)[a-e"Sb]      (2.32) 

In  order  to  obtain  the  I.C.'s  for  the  solution  of  Eq.  2.32 
make 

0°  1 

TT/.L.   n    r   -sx    /    %  j     r   _sx   -g(0)x  , 
U(t,o)  =  /  e     u(o,x)dx=/e     e  3      dx 

o  o 

=    1     n-e"(s+g(o))l  (2  33) 

s+g(o)  U  e  J  U.^J; 

which   expanded   in   McLaurins ' s    series    becomes 

„,.       ,          .          s+g(o)          (s+g(o))2  (s+g(o)  )3    +  s  +  g(o)  )4  _ 

uit,o;    -   i   -  ~y[ +  jj 47  51  ••• 

(2.34) 


or 


U(t,o)    =   hQ     (g(o))-h1    (g(l))s+h2(g(2))s2     ...  (2.35) 


00 
5    Remember  that   Lx[|£]    =      /   e"sx  |£  dx  =   2-  u(tfs! 

o 


00      .  .  .  00 


,    T       du  f      -sx  9u    ,  lim  f      -sx  8u    ,  , 

and   L     —-  =      l    e  r-  dx   =   n  /    e  ^-  dx 

x    3x  J  8x  P->oo  J  9x 

o  o 

p  p 

=   „         {e  u+s/e  u   dx} 

o  o 


=    s/    e    SX'  u   dx    -   u(t,0) 
o 

=    sU(t,s) -u(t,0) 
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Equation  2.32  for  the  derived  I.C.'s  is  known  to  have  a 
series  solution  of  the  form 

00      n  h 

U(t,s)  =  I      iiL  h  sn  =  h   -  h  s  +  =i  s  -...       (2.36) 
n=o   * 

where  h  (n=0 , 1, 2  ,  .  .  . )  depends  only  on  time. 

Substituting  Eq.  2.36  in  Eq.  2.32  and  equating  the  coef- 
ficients of  equal  powers  of  s  provides  one  set  of  ordinary 
differential  equations: 

dh 

^+  ghQ  =  f(a-b) 

dh, 

dt~+  ^hl  "  f(ho"b) 

dh 

^r2-  +  gh   =  f(nh   ,-b)  (2.37) 

dt     3  n       n-1 

The  I.C.'s  for  this  set  are  derived  setting  s=0  in  Eq.  2.36 

and  solving  as  indicated  for  the  case  of  h, (0) : 

h  (0)    u 

h,  (0)  =  ~ -H. 

1      s     s  2  22 

g     g                  g     s  +2g  s+g 
v     21    3! ;  ^   21    21   31 


(2.38) 


The  following  set  is  obtained: 


2  3 

o  o  o 

V0)   "  x  "  2T+  3T"  4T+  — 


2  3 

h    rni         X          2g°  +   39°        4g°  + 
hl(0)    =  2T  ~  3T  +   4T~  "  ST"  + 


(2.39) 


Hong  and  Larson  prove  that  h   may  be  obtained  from  the 

moment  expression   h   =  /   x   u  dx. 

n    o 
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hn  =  /   u  dx  =  xu|   -  /  x  tt~  dx 


=  u(t,l)  -  /~x  ~   dx 


The  set  of  n+1  Eqs .  2.3  7  has  n+2  unknowns  and,  for  this 
reason,  another  relationship  is  required  to  define  a  unique 
solution. 

The  approach  followed  by  the  authors  of  the  paper  is  to 
consider  the  first  moment 

1  11 

u  dx  =  xu |   -  / 
0    0 

1 

0 

1   a 
=  b  "  i  x  ^-  dx  (2.40) 

0 

Given  the  continuous  nature  of  x,  Eq.  2.40  may  be  written, 

according  to  the  mean  value  theorem  as: 

la 

hn=b-x/^dx=  b-X(b-a)  (2.41) 

0 
where 

/  xu  dx     , 

0  hl 

X  =  ~ =  p^  (2.42) 

/  u  dx       ° 
0 

The  function  b  is  obtained  from  the  above  equations  and 

it  is  given  by 

b  =  h"  -  hl  (2'43) 

Equations  2.37,  2.39  and  2.43  are  the  only  requirements 
for  an  analog  solution  of  the  problem. 

The  example  worked  out  on  the  reference  paper  considers 
only  two  equations  of  the  set  and  even  so  shows  good  results. 
As  stated  in  Section  II-B,  the  important  point  on  this 
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simulation  is  the  enormous  reduction  of  hardware  it  affords. 
This  method,  if  implemented  digitally,  leads  to  the  TSDT 
approach. 

F.   SIMULATION  BY  INFINITE  PRODUCT  EXPANSIONS  (DSTT) 

This  technique  is  based  on  the  Mittag-Lef f ler  theorem 
which  states  that  a  function  analytic  at  the  origin  and  with 
an  infinite  number  of  simple  zeros,  may  be  decomposed  in  one 
infinite  product  of  simple  factors. 

Considering  the  transcendental  function  cos  z,  it  is 
equivalent  to  the  product  (1-2z/tt)  (1+2z/tt)  (1-2z/3tt) 
(1+2z/3tt) ,  or,  if  in  Taylor's  series  form,  to  the  se- 

2     4 
z     z 
quence   1  -  ^t  +  IT  +  •   Tne  first  representation  is  char- 
acterized by  the  exact  preservation  of  the  zeros. 

Goodson  [Ref.  21]  works  out  some  examples  for  different 
types  of  p.d.e.'s.   He  considers  first  the  equations  in  the 
normal  form  and  takes  the  time  Laplace  transform  in  order  to 
obtain  an  ordinary  differential  equation  (with  s  as  parameter) 
in  state  variable  form.   Once  obtained  the  general  solution 
for  the  given  linear  boundary  conditions,  the  infinite  product 
expansion  is  applied  to  the  transcendental  functions  in  the 
solution.   The  comparison  between  the  product  expansion  and 
the  eigenvalue  or  Fourier  expansions  as  indicated  in  Figs. 
2.5a  and  2.5b,  allows  the  statements: 

(i)   Both  expansions  have  the  same  eigenvalues, 
(ii)   The  product  expansion  preserves  the  extremum  values 

of  the  exact  solution  (no  general  proof  but  it  has  been 
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Figure  2.5a.   Eigenf unction  Expansion 
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Figure  2.5b.   Infinite  Product  Expansion 

verified  in  several  cases) ,  while  the  Fourier  expan- 
sion minimizes  the  mean  square  error, 
(iii).   The  product  expansion  is  very  suitable  for  analog  sim- 
ulation due  to  the  nature  of  its  terms.   As  an  example, 
the  solution  of  the  one  dimensional  heat  transfer  equa- 
tion with  pointwise  control  at  x=0 ,  will  yield  an 
expression  as  follows: 

1 


u(i,s)  =  f(s)  n 


,  where  u(l,s)  is  the 


n=l  1+s/k 


n 


Laplace  transform  of  the  temperature  at  x=l  and  F(s) 

and  k   are  obtained  from  the  data, 
n 
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(iv)   The  poles  of  the  transfer  function  are  preserved  in  the 
case  of  the  product  expansion  and,  therefore,  the  tech- 
nique is  re commendable  for  the  simulation  of  DPS ' s 
represented  by  low-order  uni-dimensional  (otherwise  the 
analytic  computations  would  become  too  long)  p.d.e.'s. 

G.   MONTE-CARLO  SIMULATION 

Although  the  study  of  stochastic  models  is  a  little  be- 
yond the  purpose  of  this  thesis,  its  actual  importance  is 
such  that  some  considerations  are  supposedly  worthwhile. 
Two  basic  papers  were  published,  Refs.  10  and  41,  the  first 
one  describing  an  analog  computer  approach  and  the  last  one 
a  hybrid  computer  solution.   For  a  good  understanding  of 
these  papers  it  is  necessary  that  the  reader  be  familiar  with 
the  stochastic  simulation  of  systems  described  by  ordinary 
differential  equations  and  also  with  the  theory  of  Markoff 
processes.   An  adequate  theory  of  both  these  fields  may  be 
found  in  Refs.  38,  3  9  and  48. 

The  great  advantages  of  this  method  is  that  is  requires 
only  a  very  limited  amount  of  analog  hardware  or  digital  mem- 
ory and  that  it  is  only  necessary  to  make  the  computation  at 
the  points  where  the  solution  is  desired. 

H.   OTHER  SIMULATION  TECHNIQUES 

Many  other  distributed  parameter  systems  have  been  simu- 
lated in  slightly  different  ways  that  can,  however,  be 
classified  within  the  description  of  the  above  mentioned  pro- 
cedures.  Among  these  special  mention  is  made  to  Refs.  5  7 
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and  5  8  which  use  respectively  some  known  results  in  the  theory 
of  Bessel  functions  and  a  reformulation  of  the  equations  des- 
cribing a  parallel-flow  heat  exchanger  process.   By  consider- 
ing an  observer  riding  a  section  of  the  fluid,  with  the  second 
method  it  was  possible  to  reduce  significantly  the  amount  of 
hardware  required  by  the  conventional  finite  differences  an- 
alog simulation. 

I.   CONCLUSION 

After  having  described  so  large  a  variety  of  techniques, 
the  question  arises:   Which  of  them  is  the  best?, 

The  answer  is  not  a  simple  one.   Depending  on  the  avail- 
able computers  and  also  on  the  capabilities  of  each  one  (speed, 
size  of  the  memory,  number  of  components,  accuracy,  etc.), 
where  a  choice  exists,  the  most  convenient  methods  are  those 
with  which  the  user  feels  most  familiar.   It  seems,  however, 
that  sufficiently  large  analog  computers  do  not  exist,  to  per- 
mit study  of  big  or  strongly  nonlinear  systems.   In  this  case 
the  only  alternatives  are  the  hybrid  and  the  digital  computers. 

Table  I  synthesizes  the  main  characteristics  of  the  dif- 
ferent forms  of  modeling  just  described. 


46 


H 


CO 
Q 

&4 

o 

Z 
H 

a 

D 
g 
H 
CO 

fa 
O 

CO 
Q 
O 
K 
Eh 


En 
Z 


10 
•P 
P 

CD 


o 
u 


CD 

CD 
p 


CD 
5h 

cy 
x: 

4-> 

u 

rd 

cy 
p 

•H 
rH 

c 

o 
p 


p 
o 
u 

4-> 
CO 

• 

CO    CT> 

•H    C 

•H 

e  x 

<D    0) 
4J  H 

co  a 

>i-H 

CO   4-> 


CD    P 

4-1 
5h 

4-4    O 
H  4H 


T3 
CD 
CO 
P 

0) 
54 
rd 

CO 
CD 
O 
P 

cd 
5h 

<D 
IH 
4H 
•H 

I 

•H 
+J 

13 
5-1 

rci 

5-1 
O 
IH 

14-1 
•H 

Q) 

H 
X! 

rtf 
4-) 
CO 

P. 

D 


P 
•H 

54 

Q) 
rH 

rd 

O 

l 

> 

o 

p 

X! 
P 
CQ 

CD 

4-> 


CD 

t 

•H 


§ 


5-1 

rd 
iH 

P 

a. 
o 

CD 

5-4      • 
O  — 

e  p 
o 

>i-H 
<H  4-> 
0>  rd 
C  g 
•H  5-4 
CO  O 
rd  4-4 
d)  CO 
5-)  P 
U  rti 
P  V4 
H   4-1 


rd  4-4 
•H    O 
4-1 
•H    Q) 

p   CO 
•H    3 


>1 

5-4 

d  <a  «J 

5  £    CO 

O    4->    CO 
P         CD 

*    •*  o 

C    CO    0) 

p  cd  p 
> 

CJ  -H 


4->   4-> 


CO 

a 

o 
u 


CO 

C    0)  JQ 

O  "0 

•H  .— . 

4J  P  cm 

rd   co  r^ 

P   5-i 

tT'-H  . 

Q)  M-t  ft 


> 

rd 

O 
4-> 


TS 
Q) 
•H 


CO 
P 

o 

Cu-P 

cu-h 
rd  ncJ 

p 
4-4  o 

H    U 


CO 


44 
d) 


u 

> 

U 

CD 

CO 

XJ 

O 


•  CD 

co  g 

-  -H 

•  4-> 

CD 

•  5-1 
T3  CD 

•  4-) 

Cu  P 

cu 
x  e 

CD  O 

H  U 

Si 

6  P 

O  -H 
U 

4-1 
>i  CO 

5-4    O 


U 
Xi 

CT> 

•H 

X! 

CD 
£ 

4-1 
CO 

H    -rH 

Cu 

<  u 

rd 

XI 

•    £ 

0)  rd 
4->  U 
fd  TJ 
5-1 


0) 

> 

O 
4-> 

0) 

rH 

X! 

u 

•rH 


5-4  CD 
(D  Xi 
>  E-t 


CO 
P 

p  xi 

O  -H 

ax: 
5 


4-> 

rd 


4-1 

P 
Cu 


CO 

p. 

rd  4-> 

n  cu  p 
cd  x  o 

5-1    <D 
•H  4-4 

CO   4-1    O 

cu  o 

TS    P  CO 

13  CU 

CO    O  P 

■H    5-1  rH 

a  rd 
> 

0) 


c 

o 


•H   -P 
4-1  -H 


P    P  E 

rH    -H 

O  4-1  TS 

CO    P  P 

•H  rd 

0) 

X!    CD  X 

4J  XI  rd 

4J  e 
p 

CD    CO  CD 

XJ  -H  X! 

5  4-» 

rH  <-{  m 
rd  cu  a) 

e  > 

rd  54 

X  0) 


U 
■H 
4-> 


U    d)    CO 

rd 

5-1 


d) 


U 

<U 

4-> 

P 

CU 

•  6 

CO    0 

-H    O 

CO 

d>    CT> 

x:  o 

4-1    -H 

rd 

CO    C 

•H    rd 

X! 

4-J    rH 

rH 

4-1    rd 

0     & 

CO 

d) 

Cu  u 

0   0 

o 

CO    >i 

• 

5-1 

c 

Q)    0 

0 

x:  e 

•H 

4->    CD 

4-1 

6 

rd 

4-1 

4-> 

0    CD 

P 

N 

a) 

4J  -H 

e 

P    CO 

d) 

0 

rH 

rH 

CU 

CO    rH 

6 

•h  rd 

•H 

£ 

>i  CO 

4-1 

5-4 

rH 

0    CO 

p 

CD    CD 

o 

XJ   u 

•i-4 

4-1  -H 

4-1 

P 

44 

CD    G1 

■H 
Q 

££ 

fa 

H 

a 

K 

En 

fa 
O 

§ 

CO 
H 


CO 


TS  w 

CD  • 

T3  Q 

P  • 

CD  CU 

6  4H 

O  O 
U 

CD  CD 

k  a, 
>i 

Eh 


U 
■H 

rH 

o 

X) 
rd 
5-1 
rd 

04 


U 

'd  -h 

P  H 

rd    O 

XI 

5  54 

O  CD 

rH  CU 

fa  >i 


p 

rd    U 
•H 

U    rH 

•H  O 
H  X) 

O  54 
X)    CD 

rd    CU 

54     >i 

rd  B 


>H 


>i 


P 

o 

•H 
CO    CO 

p  e 

CD  CD 

6  4-1 

•H  CO 

'O  >i 

•H  CO 
P 
D 


>i 


O 
CJ 


54 

CD 

Cn 

■P 

O 

P 

^ 

Q* 

rd 

0 

S 

u 

5-4 

O 

tn 
T3  O 
•  H  H 

54    rd 


54 

O 
CT" 

T3  O 
•H  H 
54  rd 
X5    P 


54 
O  -H 

rd 

•H  -H 
54  Cn 
X)  -H 


r-{ 

rd 
4-1 
■H 

CT> 
•H 
Q 


r4 

O   rH 

rd 
T3  4-) 
•H  -H 
5h    CT> 

X3  -H 

a 


54 

°  01 

•a  o 

5h  fO 


13 
O 
X3 
4-1 


Eh 
U 
CO 
Q 


H 
Q 
CO 
CJ 


Eh 
U 
CO 
Eh 


Q 
CO 

En 


EH 

Q 

CO 
Q 


Eh 

CO 
Q 


Eh 
CO 
U 


Eh 

CO 

Eh 


fa  O 
Eh  J 

6  < 
a  u 


47 


iii.    -optimal  control  of  unidimensional  distributed  parameter 
"systems  ~~ 

a.     introduction 

In  the  sequel  of  the  general  lines  formerly  enunciated 
with  this  chapter  it  is  intended  to  create  a  physical  feeling 
for  the  optimal  control  of  distributed  parameter  systems,  by 
summarizing  different  types  of  approach  to  the  problem;  for 
simplicity  only  unidimensional  systems  will  be  considered. 
To  accomplish  this  objective  the  solutions  of  five  of  the  most 
representative  problems  are  explained  in  depth. 

As  mentioned  in  Chapter  I  the  present  work  will  not  be 
dealing  with  derivations  of  optimality  conditions  or  with  the 
mathematics  of  sensitivity,  identification,  stability,  control- 
lability and  observability,  although  making  superficial  ref- 
erences to  them.   Each  one  of  these  topics  has  in  itself  the 
potential  of  an  almost  unlimited  number  of  scientific 
dissertations . 

According  to  a  supposedly  universal  acceptance  in  this 
chapter  the  control  function  will  always  be  represented  by  u 
and  the  optimal  control  by  u*.   Similarly,  all  the  optimal 
quantities  (states  and  cost)  will  have  a  *  as  superscript. 

1.   The  Open-Loop  Control 

This  has  been  the  way  most  of  the  problems  were  worked, 
It  requires  the  knowledge  of  the  initial  conditions  and  one 
specific  field  of  application  is  in  problems  of  known 
trajectory . 
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Schematically,  the  system  representation  is  as  fol- 
lows, where  the  double  lines  represent  distributed  quantities, 
generally  in  vector  form  for  one-dimensional  systems  and  in 
matrix  form  for  multidimensional  ones.   The  optimal  control  is 
computed  and  stored  in  order  to  be  used  when  necessary. 


A    Prion 
Knowledge 

of  the  systenr 


Cost 
Computation 


7\ 


Optimization 
Procedure 


u 


:> 


Figure  3.1.   Open-Loop  Control 

Notice  that  now  the  control  is  a  function  of  the  states. 
With  exception  for  the  requirement  on  the  knowledge  of  the 
initial  conditions,  all  the  statements  in  the  last  paragraph 
also  apply  here. 

For  linear  systems  with  quadratic  cost  functions  the  opti- 
mal law  is  linear  and  an  equation  equivalent  of  the  Ricatti 
equation  for  lumped  parameter  systems  may  be  derived  and  then 
solved  in  closed  form. 

The  cost  is  a  well-defined  function  which  can  be  computed 
in  several  ways,  namely  by  numerical  methods. 

The  only  input  to  the  system  is  the  control  action,  which 
may  be  included  in  the  boundary  conditions  or  in  the  p.d.e. 
describing  it.   At  this  point  any  of  the  techniques  described 
in  Chapter  II  must  be  used. 
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The  subject  of  the  realization  of  the  optimality 
conditions  was  left  open.   Not  much  of  practical  value  has 
been  said  about  such  a  topic;  however,  a  large  amount  of 
theoretical  developments  involving  Calculus  of  Variations 
and  Functional  Analysis  is  available  in  the  literature. 

Once  again  the  approach  may  be  divided;  consider 
direct  and  indirect  methods. 

a.  Direct  Methods 

The  Direct  Methods  do  not  make  use  of  the  neces- 
sary conditions  for  optimality.   They  consider  very  often  the 
change  in  the  cost  function  (6J)  resultant  from  a  small  varia- 
tion (5u)  and,  starting  with  a  trial  solution,  they  follow 
optimization  techniques,  such  as  the  gradient  method,  until 
5J>0  for  6u>0. 

If  the  p.d.e.'s  describing  the  system  are  reduced 
to  ordinary  differential  or  difference  equation  by  techniques 
such  as  described  in  Chapter  II  (space  quantization,  time  and 
space  quantization,  eigenf unction  expansion,  transfer  function 
approximations,  etc.),  then  the  optimal  control  problem  is 
reduced  to  a  lumped  parameter,  one  for  which  there  are  well 
known  techniques. 

b.  Indirect  Methods 

The  indirect  methods  may  also  use  the  known  re- 
sults of  lumped  parameters  optimal  control.   However,  either 
using  these  concepts,  or  deriving  other  results,  they  always 
require  the  knowledge  of  the  necessary  conditions  for  opti- 
mality.  This  will  originate  equations  such  as  the  Hamilton- 
Jacobi-Bellman ,  principles  such  as  Pontriagyn's  maximum 
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principle  and  operational  procedures  such  as  the  moment  meth- 
od  and  dynamic  programming. 

Some  other  types  of  approach  are  possible,  as  for 
example  the  one  shown  by  Balakrishan  [Ref.  6]  which  includes 
constraints  in  the  form  of  a  p.d.e.  in  the  cost  function. 
2.   The  Closed-Loop  Control 

Its  general  configuration  is  shown  in  Fig.  3.2. 


U 


o 


System 
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o 
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Cost 

Computation 


t± 


Opt  imizati  on 
procedure 


Figure  3.2.   Feedback  Control 

B.   PHYSICAL  EXAMPLES  OF  OPTIMAL  CONTROL  OF  D.P.S.'s 
1.   Heating  Problems 

Suppose  a  large  heating  furnace  used  to  heat  metal 
billets  is  given  and  that  these  billets  move  uniformly  through 


This  method  was  formulated  in  general  terms  by  Kreyn 
in  1938  and  used  by  Butkowskiy.   It  is  a  good  method  within 
its  scope  of  application:   linear  systems  with  known 
eigenf unctions . 
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the  processing  zone.   The  requirement  imposed  on  the  plant, 
besides  the  cost,  is  that  the  outlet  material  must  be  at  con- 
stant temperature.   In  order  to  be  able  to  implement  a  good 
control  system  the  engineer  has  to  take  in  account  the  time  a 
given  section  of  the  material  takes  to  cross  the  furnace,  the 
temperature  distribution,  the  physical  properties  of  the 
billets,  as  well  as  its  thickness,  etc.   Intuitively,  the  pos- 
sible control  laws  would  take  the  form  of  changing  the  speed, 
or  the  temperature,  or  both. 

There  are  numerous  constraints  to  be  taken  in  considera- 
tion, as  for  instance  in  the  furnace  temperature  T,(x,t) ,  such 
that  C~<T,     (x,t)  <  C,,  where  C,  and  C«  are  constants  and  x  is 
the  distance  from  origin.   Also  the  interaction  between  adja- 
cent zones  of  the  furnace  jnust  obey  the  inequality 

8T1(x,t) 

— 5 <C0  ,  the  billets  temperature  T_(x,y,t)  must  be  less 

dX  j  Z 

than  a  constant  (here  y  is  measured  in  the  direction  of  the 

dT2(x,t) 
thickness)  and,  similarly,  * <C„  . 

2 .   Heat  and  Mass  Transfer  Problems 

One  example  of  this  type  of  problems  is  the  process  of 
drying  of  a  moist  material  in  a  continuous  drying  unit.   In 
this  case  the  control  action  consists  in  compensating  the 
changes  in  the  composition  of  the  moisture,  porosity  of  the 
material,  layer  thickness,  velocity  of  the  material  inside 

the  heater,  etc.   The  constraints  may  assume  different  as- 

'  i  8u  i 

pects,  a  typical  example  being  m<u(x,t)<M,  |— |<C,  and  the 

temperature  inside  the  material  also  less  than  a  constant. 
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3 .  Separation  of  Liquids  and  Gaseous  Mixtures 

When  using  big  reactors  their  distributed  nature  must 
be  considered  and  once  again  the  D.P.S.'s  techniques  are  re- 
quired.  The  problem  may  be  formulated  in  three  different 
ways,  in  terms  of  the  control  action: 

(i)   Control  of  the  source  heat  to  obtain  the  highest  con- 
centration of  the  substance,  assuming  a  constant  output 
rate  ; 
(ii)   Control  of  the  source  heat  to  obtain  a  maximum  product 
output  for  a  given  concentration  of  the  mixture;  and 
(iii)   Minimization  of  the  control  costs  for  a  given  concen- 
tration and  output  rate  of  the  product. 

4 .  Automatic  Control  of  Large  Hydroturbines 

In  this  case,  the  distributed  parameter  nature  of  the 
pipeline  supplying  water  to  the  turbine  has  to  be  considered. 
The  control  of  the  turbine  rotor  is  realized  with  a  valve, 
but  while  this  action  should  be  fast,  in  order  to  minimize 
abrupt  changes  in  velocity  when  the  load  is  removed,  it  has 
to  be  limited  by  the  resulting  changes  in  pressure  on  the  pipe- 
line walls,  which  cannot  exceed  the  safety  value. 

5 .  Gas  Transfer  Through  Long  Pipelines 

The  problem  consists  in  locating  and  controlling  com- 
pressors along  the  pipeline  such  that  the  variations  in  the 
output  pressure  are  minimized.   Considering  only  one  compres- 
sor situated  at  the  origin  (single  pointwise  control)  the 
boundary  conditions  will  be  of  the  type 

p(o,t)  =  u(t) 

vU,t)  =  v  (t)  , 

a 
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where  p(o,t)  is  the  pressure  at  x=o  and  v(£,t)  is  the  gas 

velocity  at  x=£.   This  velocity  equals  v  (t)  which  in  turn  is 

a 

a    function    of   the    load. 

The    functional    to  be   extremized   is    of   the    type 

T 
J  =    !    |p*(t)     -   p(£,t)Y|dt    ,    y*1 
o 

The  control  of  systems  such  as  those  just  described 
is  by  itself  a  challenging  task  due  to  the  complexity  of  the 
equations  involved.   Unfortunately,  this  task  is  made  even 
more  difficult  by  the  fact  that  some  of  the  parameters  of  the 
partial  differential  equation  are  not  accessible,  which  im- 
plies the  need  for  the  use  of  parameter  identification 
techniques . 

C.   FIVE  REPRESENTATIVE  EXAMPLES 

A  brief  analysis  of  the  contents  of  this  section  is  given: 
Example  No.  1  [Ref.  5]  is  an  open -loop  indirect  method  appli- 
cable to  a  large  variety  of  D.P.S.'s  for  which  a  pointwise 
control  is  desired.   It  is  based  on  a  space  discretization 
technique  that  leads  to  a  set  of  ordinary  differential  equa- 
tions to  which  the  methods  of  optimal  control  of  lumped  para- 
meter systems  can  be  applied. 

Example  No.  2  [Ref.  5]  has  similar  characteristics  to  the  pre- 
ceding example  but  is  a  multiple  pointwise  control  problem. 

Example  No.  3  [Ref.  56]  is  the  last  of  the  open-loop  problems 
presented  here.   It  is  also  an  indirect  method  and  may  be  used 
whenever  the  partial  differential  equation  may  be  reduced  to 
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the    form  Q (x, s) =G(x,s) U (s) ,    where   Q(x,s)    and   U(s)     are    respec- 
tively  the    transforms    of   the    output    and   of    the    input. 

Example   No.     4     [Ref.    29]    is    one    of   the    few    closed-loop   methods 
available    in   the    literature.       It    applies    to  parabolic   differen- 
tial equations   with   quadratic    cost-function.      The    control    is 
multiple   pointwise. 

Example   No.    5     [Ref.    28]    presents    again    a    closed-loop   technique. 
It   applies    to   partial    differential  equations   with  known    anal- 
ytical  solution    and  when   the    functional    cost    is    quadratic.       In 
the   example    shown    there    is    only    a   single    control  but    a  multiple 
pointwise    control   is    also  possible. 

1.       Open-Loop    Minimum  Time    Optimal    Control    of    a   Heat 
Transfer   System   Using   Space    Discretization    and 
Pontryagin's    Maximum   Principle 

Consider   the   heat   equation 

2 

!§■  =   k|-#   ,    0<x^L    and    0^t<T  (3.1) 

9t  9x2 

The   boundary    conditions    are 

-    4^|x=Q    =   b[u(t)     -   q(0#t)]  (3.2) 

where    a    and   b    are    given    constants,    and 

■^2.1  =    0  (3    3) 

9x'x=L         U  U#jj 

The  initial  condition  is 

q(x,0)  =  qI>c> (x)  (3.4) 

The  pointwise  control  is  bounded:   |u(t) | <M 
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One  way  of  reformulating  a  minimum  time  problem  is 

to  choose  a  convenient  value  for  the  final  (tf)  such  that  at 

this  instant  the  deviation  from  the  optimal  trajectory  is 

minimum.   An  adequate  cost  function  is  given  by 

L 
J  =  /  |q*(x)  -  q(x,t)|Ydx  ,  y>l  (3.5) 

o 

Discretizing  equations  3.1  through  3.4,  the  following 

expressions  result  for  L=NZ : 


q2-q1      q1-g0 

ql   =   k.     I  1/2 


dt  I 

qitl"qi        qi"qi-l 


dq, 


=  k — 5 ,    i=2,3,...N-l         (3.1a) 


dt  I 

qN+l_qN  qN"qN-l 


dq 


n   =   .        1/2 l_ 


dt  I 


ql"q0 


-    a      z>2      =   b(u-qQ)  (3.2a) 

This    is    clearly   illustrated   in    Fig.    3.3. 

The    above    discrete   equations    can   be   written    in    dimension- 
less    form: 

%  -  m  u  +  m  qi  '    for  8  =  !r  <3-2b) 
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(N-2)l  (N-DI   L  x 


)  . 
Figure  3.3.   Method  of  the  Straight  Lines 


dq 


=  2qn  -  3q,  +  q9  ,   T  = 


'0 


•1    ^2 


2t 

12 


dx 

dq± 

d~  =  qi_X  "  2q±  +  q.+1  f   i=2,3,..-,N-l 

dq 
__n_ 

dT  -  %_i   qN 

Substituting  (3.2b)  in  (3.1b)  the  result  is 

dql    23      2+33 

dT  ~  2+3U    2+3   ql   q2 


dq 


d 

dq 


7^  =  qi-l  "  2qi  +  qi+l  '   i=2'3'*"  'N_1 


n 


dx     qN-l    qN 
which  is  easily  implemented  on  an  analog  computer. 


(3.1b) 


(3.1c) 
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The  given  p.d.e.  may  also  be  simulated  easily- 
remembering  that 

3t    re  8X2 


(3.6) 


gives    the    voltage    distribution    in    a  non-inductive    transmission 
line   that   has    r   and    c   as    resistance    and    capacitance   per   unit 
length.      Therefore    the    total    resistance    and    capacitance    are 
R=rL   and   C=cL. 


R- 


Qu(t) 


Figure  3.4.   An  Analogic  Simulation  of  a  p.d.e. 


The  equivalent  m  sections  circuit  is  represented  in 

R  I? 

Fig.  3.4  where  g-  =  2£N  ,   rc  =  =- 

K-j  a 

and  the  load  R  =  °°  because  of  Eq.  3.3. 


(3.7) 


Returning  to  Eq.  3.1c,  suppose  that  the  following 
initial  conditions  are  given 
q0(0) 


(3.8) 


becomes 


q^O)  =   qQ(0)  (^-^l)     ,   i=l,2,...,N 
Equation  3.5  reformulated  for  the  discrete  case 

n 


J  =   I    q  Z(t  ) 
i=l  x        r 
which  is  to  be  maximized. 


(3.9) 
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The  system  (3.1c)  is  representable  in  the  linear  form 


as 

x(t)  =  A(t)  x(t)  +  Bu  (3.10) 

and,    as    a    consequence,    the   necessary    condition    for   optimality 

from  Pontryagin's    principle   becomes    also  the    sufficient 
7 


condition. 


Writing   the   Hamiltonian 


8 


t  ,23  2+33      v 

H   =    -    1    +    p1(2--u-rFrq1)     +    q2 

n-1 
+    J2Pi  Cqi-l~2qi+qi+l)  +PN  (%-l-%J 


(3.11) 


where  the  p's  are  the  Lagrange  multipliers. 

As  it  is  easily  seen,  given  that  [  u  |  <M  ,  the  Hamil- 
tonian has  a  maximum  for 

(3.12) 
Looking  now  at  A  of  Eq.  3.10: 
2+36 


u(t)    =   M   sign    P-,  (t) 


A      = 


2  +  3 

1 


10         0 0         0 

-2         10 0         0 


(3.13) 


0  0         0         0  -2         1 

0  0         0         0  1-1 

It  should  be  noticed  that  it  is  symmetric  and  therefore  its 
eigenvalues  are  real.  It  also  can  be  shown  that  the  eigen- 
values   are    all   negative.      This    last    result   shows    that    the 


Reference    65,    p.    2  34. 


Reference    37,    p.    188. 
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system  moves    to   the    origin    as    t-*°°   and   this    implies    the   exis- 
tence   of   a   unique    optimal    control.      Because    the    eigenvalues 
are    real    and  negative    the    optimal    control   will    change    its 
value   n   times,    as    illustrated   in   Fig.    3.5    for  n=4 .      The    ini- 
tial    and    final    conditions    are   q. (0)=-l    and  q. (T)=0    for 
i=l,2,3,4, 


Figure    3.5.      The    Optimal    Control    Law    and   the    Correspondent 
Output. 


2 .       A   Simple   Example   with    Open-Loop    Multiple    Pointwise 
Control 

Consider   the   heat   transfer   equation 

a(x,t)     |2.  +    a(x,t)v(t)    |SL  +   q   =    u(x,t)  (3.14) 

between   a   medium   at   temperature    u    and    a  body   moving  with   veloc- 
ity  v >.0 .      The   heat   takes    place   between    the   x   origin    and   x=L. 
The    remaining    information    is:       a   and   v   are   known    and 

(3.15) 
(3.16) 


q(x,0)    =   qI    c    (x) 


q(0,t)     =    0    ,  0<t<T 
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The  problem  consists  in  finding  the  control  u(x,t) 

such  that  the  functional 

T 
J  =  /  |q*(t)  -  q(L,t)  |Ydt,  Y>1  (3.17) 

0 

is  minimized  or  equivalently ,  such  that  the  time  integral  of 
the  temperature  deviation  of  the  output  material  from  a  given 
law  is  minimized.   The  following  constraints  are  imposed: 

<Z      <  u(x,t)     <   C, 

au  (3*18) 

4  8x  3 

Discretizing   the    system   in   n    linear  elements,    Eq. 
3.14   becomes 

dq .  q . -q . _ , 

a.  (t)    -htF-  +    a.  (t)v(t)       \    1        +  q.    =   u.  (t)  (3.19) 

i  dt  l  I  ^i  l 

But,  by  Eq.  3.16,  q  (t)  =  0  and  therefore  it  follows 

dq .  u. 

dT~   =  Bqi-1  +  aiqi  +  a^"  '   i=1'2" ■ •  'N  (3.20) 

i 

where 

a  V(t)      ,  1        V(t)  ,->   -  -,  x 

3  =  ___  and  a.  =  -  ___  -  __  (3.21) 

In  order  to  have  the  term  with  u.  free  of  other  time  functions 

l 

■ 

define  q.  =  a.q.  and  rewrite  Eq.  3.21  as 


dq!. 

dt     ^1-1     i^i     l 

The  discrete  form  of  the  constraints  becomes 

C0  <  u.  (t)  <  C, 
2  —  l     —   1 

C,  <    u.,,  -  u.  <  C- 
4  —   l+l     l  —   3 


(3.22) 


(3.23) 
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and  the  cost  takes  the  form 

T 
J   =   /  |q'*  -  q'(t)|Ydt,   Y>.1  (3.24) 

0  n 

The  corresponding  Hamiltonian  is 

N 
H(p,q,u)  =  P0|q,*(t)-q^(t)|Y+<X  p±  (Bqll^+cuq  |+u± )     (3.25) 

and,  as  usual,  the  costate  equations  are  obtained  by 


dt  9q' 


or 


dp* 
dt 

=      0 

dp* 
dt 

=    - 

dp* 

^n 

aipi  "  Spi+1  '   i=l/2,««-,N  (3.26) 

dt    "  "  ^?olq'*  "  qN|Y_1  sign(q'*-q^)  -  a^ 

As  it  may  be  immediately  seen,  for  p.  =  const.  H 

takes  the  maximum  value  for  a  corresponding  maximum  of 

N 
F(p.,u.)  =   £  p*  u.   satisfying  the  constraints  in  Eq.  3.23. 
i=l 

If  L=2£  it  is  possible  to  draw  one  elucidative  two- 
dimensional  picture  showing  the  effect  of  the  constraints  on 
the  optimal  control.   Equations  3.26  together  with 

4rj—  =   ?=-  ,  the  control  constraints  and  H(q*(T),  u*  (T)  ,  p*(T), 

d  t         dp  ~         ~         ~ 

*"  q 
T) 6T  =  0     lead  to  the  complete  solution  of  the  problem. 


9 

Reference  37,  p.  2  33. 
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d 

U2 

0, 

C2 

c, 

u, 

/ 
/ 
/ 

/ 

Figure  3.6.   Geometrical  Interpretation  of  the  Control 
Constraints. 


Butkovskiy  [Ref.  5]  describes  two  other  methods  used 
to  find  the  switching  times  of  the  control  in  minimum  time 
problems.   They  are  the  method  of  harmonics  and  the  method  of 
parabolic  approximations.   Both  give  better  accuracy  than  the 
straight  lines  approximation  illustrated  in  Fig.  3.3  and  they 
are  mentioned  here  as  a  reference  to  the  interested  reader. 

Finally,  since  there  are  few  references  to  multi- 
dimensional systems,  it  seems  useful  to  refer  to  the  above 
mentioned  work  for  an  explanation  of  a  simple  way  of  reducing 
the  problem  of  optimal  heat  of  a  sphere  or  of  a  cylinder  to 
the  one -dimensional  problem  using  the  method  of  harmonics. 
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3.       Sakawa's    Solution    of    One    Open-Loop    Optimal    Heat 
Transfer   Problem 

Sakawa's    solution    [Ref.    56]    is    one   of   the    most    com- 
prehensive   and    complete    papers    in    the    field   of   practical    ap- 
plications   of    distributed   parameter   systems. 

Given   the   normalized  parabolic   equation 

2 

9    q(x,t)  3q(x,t)  ,  . 

„    2  3t  K6'z  n 

oX 

for  0<x<l  and  0<t<T  with  initial  and  boundary  conditions 
q(x/0)  =  0 

^Bx'^  'x=0  =  a^q(°'t)  -  v(t)}       .        (3.28) 

3q(x,t) I     _  0 
3x   'x=l 

the  control  u(t),  0<t<T,  is  a  fuel  flow  which  controls  v(t), 

the  temperature  of  the  medium. 

Equation 

Y  ^t^"  +  V(t)  =  U(t)  (3.29) 

expresses  the  relation  between  u(t)  and  v(t)  and  is  physically 
equivalent  to  a  first-order  lag. 

The  objective  is  to  minimize  the  functional 

1  2 

J(u(t))  =  /  (q*(x)  -q(x/T)}^dx  (3.30) 

0 
where  q*(x)  represents  the  desired  temperature  distribution 

at  t=T. 

Applying  Laplace  transform  techniques  to  Eqs .  3.2  7 

and  3.29  and  solving  for  the  given  boundary  conditions, 

Q(x/s)  may  be  written  as  Q (x, s) =G(x, s) U (s) . 
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From  this  equation,  by  the  convolution  theorem 

t 
q(x,t)  =  /  g(x,T)u(t-T)dT 
0 

(3.31) 
t 
=  /  g (x, t-x) u (x)  dx 
0 

Note:   All  the  procedures  that  follow  may  be  applied  not  only 
to  the  heat  equations,  but  to  any  partial  differential  equa- 
tion reducible  to  the  form  of  Eq.  3.31. 

N  (x  s) 
Rewriting  G(x,s)  as  the  quotient   G(xfs)  =   '  ' 

it  is  easy  to  use  the  theorem  of  the  residues  and  obtain 

°°   N  ( x ,  s  .  )   s  •  t 

g(x,t)  =  I —  e  X  (3.32) 

i=0  M1 (s. ) 

where  the  m'(s.)  is  the  derivative  of  M(s)  at  s=s. 
1  1 

By  the  use  of  standard  minimization  techniques  des- 
cribed in  detail  in  the  reference  paper,  from  Eq.  3.30,  the 
optimality  condition  can  be  derived: 

1 
/  (q*(x)  -  q(x,t) }g(x,T-x)dx  =  0  (3.33) 

0 

1 
Defining  f(x)  =  /  q* (x) g (x,T-x) dx  (3.34) 

0 

it  comes  out,  by  bottom  Eqs .  3.31  and  3.33: 

T       1 
f(x)  =  /  u(y)  /  g(x,T-x)g(x,T-y)dxdu  (0<x<T)        (3.35) 
0       0 

where  u  is  a  time  dummy  variable. 

The  integral  on  the  right 

1 
yU,u)  =  /  g  (x,T-x)  g(x,T-u  )  dxdu  (3.36) 

0 
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is    a   symmetric  kernel        and   as    such   has    some    important 

properties : 

(i)   There  exists  at  least  one  eigenvalue,  X .  ^  0  3 

T 
<f>    (x)    =  X       J   y(x,y  )  cfr.  (y  )  du    where    <f> .  (x)  (0<x<T)    is    the 
0  1 

eigenf unction    corresponding   to  X  .     . 

(ii)       The   eigenf unctions    are    mutually   orthogonal,    i.e., 

T 

/    <j>.  (x)  <J>.  (x)dx    =    0    ,       for  X  .    i-  X  .    .' 
0      i  D  13 

(iii)   The  necessary  and  sufficient  conditions  for  existence 

00 
of  solutions  of  Eq.  3.35  is  that   >  X.  c.  must  be  con- 

i=l    X      X 
T 


vergent,    where    C.    =    /    f  (x)  <J> .  (x)  dx    .      The    optimal 

1        0  L 

control  will   be    given   by: 


u*(t)    =      I      X    c.  <J>.  (x)     .  (3.37) 

i=l      -1-    x    x 

From  Eq.  3.37  it  can  be  seen  that  in  order  to  obtain  the 
optimal  control  it  is  necessary  to  solve  also  Eq.  3.35  which, 
in  general,  is  a  hard  task.   Also  it  happens  that  sometimes 
the  solution  of  Eq .  3.35  does  not  account  for  control  restric- 
tions.  For  these  reasons  Sakawa  develops  a  numerical  integra- 

1 

tion   of   Eqs.    3.31    and    3.32    obtaining 

J(u)     -       I    C.(q*    -      I    a         u    )2  (3.38) 

i=0  j=0      J      J 


See    Ref.    11,    p.    115,    for   detailed   definition    and 


properties . 
11 


3  means  "such  that." 
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where  the  C.'s  and  a. ,'s  are  obtained  in  a  straightforward 
manner.   A  way  of  minimizing  Eq.  3.38  is  by  quadratic  program- 
ming/ for  which  Sakawa  gives  appropriate  references.   However, 
because  a  linear  programming  technique  is  much  easier  to 
implement,  he  considers  another  cost  function: 

1 
J(u(t))  =/  |q*(x)-q(x,t) |dx  (3.39) 

0 

From  this  he  derives,  in  a  similar  fashion  as  for  Eq.  3.38, 

the  following  equation: 

S    ,  *     9 
J(u)  -  I    C. |q.  -  I    a   u  |  .       (3.40) 

i=0  L   x    j=0  1J  : 
The  control  constraints  are  taken  in  account  in  the  linear 
(or  nonlinear)  programming. 

In  the  sequence  of  this  last  procedure  the  author  ob- 
tained the  following  curves  for  the  parameters 

a  =  10,    Y  =  0.04,    q*(x)  =  .2,    n  ■=  20 

Summarizing,  Sakawa' s  method  is  considered  to  be  quite 
valuable  for  single  pointwise  control  of  unidimensional  sys- 
tems.  Although  it  is  not  impossible  to  generalize  this  meth- 
od for  multiple  control  using  multiple  Laplace  transforms  and 
the  superposition  principle,  this  seems  to  be  a  huge  task  and 
therefore  the  method  is  not  recommended  for  such  cases. 

Quite  recently,  Chang  [Ref.  9,  August  19  70]  developed 
one  algorithm  based  on  a  modified  steepest  descent  method  and 
solves  the  same  problem  as  above  obtaining  very  good  agreement 
with  only  two  iterations. 
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Figure    3.7.       Optimal   Control 
(T   =    .2)  . 


Figure    3.8.       I.      Temperature 
Distribution    for   u(t)    =    1. 

II.      Temperature    Distribution 
for   u* (t) . 
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Figure  3.9.   Optimal  Control 
(T  =  .4)  . 


Figure  3.10.   I.  Temperature 

Distribution  for  u(t)  =  1. 

II.   Temperature  Distribution 
for  u*(t) . 


68 


4 .   Greenberg's  Solution  of  a  Closed-Loop  Multiple 

Pointwise  Optimal  Control  Problem  with  Quadratic 
Cost  Function 

a.   Pointwise  Optimal  Control 

The  mathematical  derivations  used  by  Greenberg  in 
his  paper  [Ref.29]  are  hard  to  follow  if  the  reader  does  not 
have  a  deep  understanding  of  Functional  Analysis  and  theory 
of  operators.   However,  the  consequences  of  this  method  are  so 
important  and  its  implementation  is  so  simple,  if  compared 
with  the  complexity  of  the  derivations,  that  it  is  worthwhile 
to  analyze  its  contents  in  some  detail. 

The  typical  approach  to  the  problem  represents  a 
distributed  control  law  by  a  finite  number  of  lumped  controls, 
but  this  is  not  applicable  to  a  feedback  control.   Also  the 
modal  control,  by  the  reasons  already  pointed  out  in  Chapter 
I,  is  difficult  to  realize  if  the  eigenvalues  are  hard  to  ob- 
tain or  if  they  converge  slowly.   As  it  will  be  shown,  the 
present  method  does  not  suffer  such  drawbacks. 

The  author  starts  with  the  knowledge  previously 
demonstrated  that  systems  described  by  parabolic  differential 
equations  with  quadratic  cost  functions  have  an  optimal  feed- 
back control. 

Writing  the  parabolic  differential  equation  in 

the  form: 

q(t)  =  Aq(t)  +  B(t)u(t)  ,  q(o)  =  q   ,  (3.41) 

o 

where   A   is    a  partial    differential    operator    (generally   of   the 

type    ■s — j  )     and   B  (t)    is    a   time   variant    control    operator   defined 
x 

for    D   <    t    <   T,    the    functional    cost    is    defined    as 

69 


T  12 

J  =  /  [<Q(t)q(t)  ,q(t)  >+uT(t)R(t)u(t)  ]  dt 
0  ~ 


(3.42) 


R(t)  is  a  kxk  positive  definite  matrix  and  Q ( t)  a  positive 
operator;  k  represents  the  number  of  desired  control  points. 
After  non-trivial  manipulations  the  optimal  con- 
trol law  is  found: 


u*(t)  =  -  R~1(t)B(t)/  k(t,c)q(t,C)d£ 

D 


(3.43) 


where  k.(t,x)  =  K(t,x,x.) ,  i  =  1,2, k  with  K(t,x,£)  the  solu- 
tion of  a  given  integro  -  differential  equation,  and 
B(t)  a  kxk  diagonal  matrix,  with  easily  computable  elements. 

The  block  diagram  correspondent  to  the  above  men- 
tioned description  is  represented  in  Fig.  3.11. 


u(t,x)-^(t)u   Dynamical  System 
^  =  Aq  +  u(t.x) 


Pointwise 

Control 

Operator 


<^ 


-R(t)B(t) 


£ 


q(t,x) 


k(t,xX)dx 


k    vector 

distributed 
quantity 


Figure    3.11.      Block   Diagram    for   Greenberg's    Method 


12 

The    symbol    <x,y>   means    inner  product    of   the    vectors 

x   and  y . 
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B  (t)  is  defined  using  spectral  decomposition  of  operators, 

a  subject  that  will  be  treated  in  detail  in  the  next  chapter. 

b.   Infinite  Terminal  Time  Problem 

The  above  description  can  be  quite  simplified  in 

the  case  of  a  terminal  time  problem  because  the  operators  B 

o 

Q,    R   and  K   the    feedback    operator   become    time    invariant. 
Represent 

Q(x,c)    =   yT(x)9Y(C)  (3.44) 

where   V(x)    is    the   vector   of   the    first  N_  eigenfunction    of   A 
and  Q    a  N   x  N   positive    definite   matrix.       Define   K    as    the  posi- 
tive  definite    solution    of    the  matrix    Riccati   equation 
-AK-KA+KVTB    R-1    B    V  K    -    Q    =    0,    with    A    a   N    x   N    diag- 
onal   matrix   such   that    A. .    equals    X.    of   A.       Also   define    V   as 

a  k   x  N   matrix  whose    ij         element    is    V.  .    =    v. (x.) .      The    optimal 

J  ij     j      l         r 

control  is  given  by 

u*(t)  =  -  R"1  B  V  K  /v(C)q(t,C)dC 

D  (3.45) 

=  R_1  B  V  K  q(t) 

where  q(t)  is  the  N  vector  of  modal  coefficients  of  the  state 
distribution  q(t,x)  . 

Due  to  Eq.  3.44,  if  the  optimal  control  exists  it 
cannot  be  improved  by  feeding  back  more  than  N_  modes.   An- 
other important  conclusion,  is  that  whenever  (QN+i~Q^  i-s    a 
positive  operator  it  follows  that  K^+,-K>0.   This  can  be 
stated  as  "monotone  approximation  of  the  state  weighting 
operator  originates  monotone  approximation  of  the  optimal  feed- 
back operator." 
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Greenberg  finishes  his  paper  with  a  very  simple  and 

illustrative  example. 

5 .   Graham's  Solution  of  a  Close-Loop  Pointwise  Optimal 
Control  Problem 

J.  W.  Graham's  paper  [Ref.  28]  presents  two  different 
approaches  to  the  solution  of  an  optimal  control  problem  char- 
acterized by  a  given  parabolic  differential  equation  with  a 
frequent  type  of  initial  and  boundary  conditions.   The  author 
uses  two  types  of  approach  to  the  problem;  the  first  one  deals 
with  the  matrix  Ricatti's  equation  and  the  second  one  con- 
siders Kalman's  equation,  together  with  root-locus  techniques, 
constituting,  in  a  certain  way,  like  a  smooth  transition  be- 
tween the  classical  and  the  optimal  control  problems. 

Some  details  in  the  derivations  are  very  hard  to  fol- 
low due  to  gaps  in  the  explanatory  theory;  for  this  reason 
additional  information  will  be  given  in  an  attempt  to  clarify 
the  points  where  more  omissions  were  found.  The  nomenclature 
used  is  consistent  with  the  remaining  of  this  chapter,  but  the 
reverse  of  the  one  in  the  paper  regarding  the  control  and  the 
temperature . 

The  p.d.e. 

2 

|§.(x,t)  =  v2  |-4  ,      t>0  ,    0<x<a  (3.46) 

3t  3x2 

is  given  with  initial  conditions 

q(x,0)  =  QQ(x)  (3.47) 

and  boundary  conditions  (B.C.'s) 

H  (0't}  =  ° 

K|3.(a,t)  =  f(t)  -  q(a,t)  =  u(t)  (3.48) 
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where  K  is  the  thermal  conductivity  and  f(t)  is  an  external 
available  temperature  input.  All  the  other  symbols  are  de- 
fined as  before. 

The  cost  function  is  defined  as 

oo 

J  =  /  [q(a,t)-qss(a)  ]2  +  R[u  (t)  -u^]  2dt  (3.49) 

where  the  subscript  ss  represents  steady-state  and  R  is  a 
positive  constant.   The  physical  meaning  of  the  equation  is 
obvious . 

The  analytic  solution  of  the  p.d.e.  for  the  given 
boundary  conditions  and  piecewise  constant  control  (u(kt) )  is 
given  (although  not  necessary) : 

q(x,t)    T,  rv  t  ,  3x  -a     2    r  (-1)      ,   v  n  tt  t,    n-nx 

ulkrT    =  Ka{"lT  +  — —  "  ~2     ^~T-  exp[ 2 ]cos-T- 

6a      7T   n=l  n  a 

(3.50) 
for  [kT<t<(k+l)T]  . 

After  this  the  author  states:   "The  original  sampled- 
data  equations  in  discrete  form  are  transformed  to  the  normal 
or  diagonal  form  in  a  straightforward  manner  and  from  these 
equations  a  continuous -time  set  of  equations  is  readily  ob- 
tained.  Hence  the  system  equations  may  be  written  as  the 
following: 

W(t)  =  AW(t)  +  B  u(t)  (3.51) 

Q(t)  =  C  W(t)  (3.52) 

where  W(t)  is  an  N  dimensional  column  vector,  A  a  NxN  diag- 
onal matrix,  B  and  N-dimensional  column  vector,  Q(t)  the 
temperature  of  the  system  at  x=a  and  C  is  an  N-dimensional 

vector. " 
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The  way  Eqs .  3.51  and  3.52  were  obtained  does  not 
seem  so  obvious  as  the  author  pretends.   One  possible  tech- 
nique that  can  be  used  to  find  them  is  the  Bubnov-Galerkin 
method  described  in  Appendix  A. 

Back  to  the  original  problem,  the  functional  J  can  be 
further  simplified  by  defining 

Q*(t)  =  Q(t)  -  Q, 


ss 


u*(t)  =  u(t)  -  u 


ss 


from  which 


W*(t)   =  AW* (t)  +  B  u*(t) 
QCt)    =   C  W*(t) 
with  initial  conditions 

Q*(o)  =  -Q, 


'ss 


W*(o)  =  -W 


ss 


and   it    comes    out 


J   =      /    [Q*(t)2    +    Ru*(t)2]dt 
0 


(3.53) 


(3.54) 


(3.55) 


(3.56) 


or 


J   =    /     [W*(t)TSW*(t)     +    Ru*(t)     ]dt 


(3.57) 


where 


S      = 


C    C 


Vl 


c  c 

U1U2 


c  c 


'l0® 


c2°n 


'N 


(3.58) 


and   the   Ci's    (1=1,2, N)    are    the   elements    of   the    C   vector. 

Note   that   in   the    above   equations    the    superscript    (*)    is    a 
mean    for    identification    of    some    variables    after 
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transformations  on  them,  rather  than  being  the  usual  con- 
vention for  optimal  function. 

At  this  point  the  author  considers  two  different  ap- 
proaches to  the  computation  of  the  optimal  control. 

a.   Matrix  Riccati  Equation  Method  of  Solution 

Because  of  the  linearity  of  the  given  system  and 
also  because  a  quadratic  cost  function  was  chosen,  it  is 
known  that  the  computation  for  the  optimal  control  may  be  re- 
duced to  the  solution  for  P  of  the  Riccati  equation  [Ref.  3, 
p.  752-775] 

P(t)B  BTP(t) 
P(t)  +  S  - "V +  P(t)A  +  AP(t)  =  0        (3.59) 

with  boundary  condition  P(tf)  =  0.   This  boundary  condition 

comes  from  the  definition  of  J.   If  J  had  another  quadratic 

T 
term  outside  the  integral,  i.e.  the  term  W*  GW* ,  the  boundary 

condition  would  be  P(tf)  =  G. 

The  optimal  control  law  is  given  by 

u*(t)   .  =  -  R-1BTPW(t)  (3.60) 

opt      „      „   

where 

P  =  f-im  p(t)  (3.61) 

t  f -*-°°   ~ 

whenever  the  system  is  controllable  and  has  P(tf  =  0). 
P  is  obtained  as  the  solution  of 

-  PA  -  AP  +  PBR~1BTP  -  S  =  0,  (3.62) 

or  by  solving  Eq.  3.59  backwards  in  time  from  the  known 
condition  P(°°)  =  0  until  convergence  is  reached.   This  last 
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method  can  be  implemented  very  fast  in  a  digital  computer 
using  the  algorithm 

P(t+A)*P(t)+A{-P(t)A  -  AP(t)+P(t)BR_1BTP(t)-S}  (3.63) 

The  graphical  interpretation  of  the  procedure 
is  given  in  Fig.  3.12 


0 


P(t)=P 


Steady  State- 


-u 


PCtf  )=0 


T— Tpo      t 
Transient  r — 


Figure  3.12.   Graphical  Interpretation  of  the  Solution  of 
Equation  3.67. 


It  is  important  to  note  that  P(t)  is  independent 
of  the  states  and  can  be  computed  once  J  is  specified,  before 
the  optimal  system  is  started. 

Once  u* (t)     is  obtained  from  Eq.  3.60,  the  final 

opt  ^ 

result  u(t)     is  given  by  u*(t)     +  u 

opt  opt     ss 

b.   Kalman's  Equation  Method  of  Solution 

The  foregoing  method  could  be  applied  without  the 
need  for  reducing  the  state  equations  to  the  diagonal  repre- 
sentation.  This  is  not  the  case  in  the  present  paragraph, 
because  the  nice  results  derived  from  Kalman's  equation,  and 
explained  in  deep  detail  in  Chapters  VII  and  VIII  of  Ref.  60, 
are  much  more  easily  used  for  the  diagonal  form. 
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Figure  3.13.   The  Global  System 

Suppose  an  arbitrary  i  row  of  an  n  dimensional 
state  equation: 


w.  =  X . w.  +  b . u  . 
1     11     l 


(3.64) 


Taking  the  Laplace  transform  and  ignoring  the 
initial  conditions  for  the  sake  of  clarity,  the  distributed 
system  can  be  represented  as 


b, 
s-X, 

c, 

b 
S"X2 

c2 

-4- 

U(s) 

bi 

L            A 



i 
i 
i 
■ 

1 

1 

I 

W), 

Vl 

1 
1 

J 

bN 

S~Xk 

u 

Figure  3.14.    Approximate  D.P.S.  Transfer  Function 

(becomes  exact  as  n^-°°) 
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GN-1  cN-2 

CM  Q*(s)  ^-l3  +    aN-2S  +    —    +    aQ 

^s;         U*(s)    "        (s-L)(s-L)— (s4J 


(3.65) 


But,    as    shown    in    Ref.    60,    from  Eq.    3.6  6    also 
known    as   Kalman's   equation    it   can   be    derived   that 


|1  +    G(s)H       (s)  |2    =    l+I^M-l2 
e<3  /R- 


where 


(3.66) 


|1+G(S)H       (s)T    =    (l+G(s)H       (s))(l+G(-s)H       (-s)) 

|G(s) |2  =  G(s)  G(-s) 
and  that  the  optimal  control  is  givey  by 


(3.67) 


U*      L(s)    =   r-hW*(s)  (r=0    in    the    present    case)  (3.68) 

opt 


such   that 


H       (s) 

eq 


hW*(s)  hxw*(s)    +   h2W*(s)    +    - 


-    +    Vn(s) 


3.15b. 


Q*(s) 


c.wJCs)    +    c2w*(s)    +    - 


-    +    CNWN(S) 


(3.69) 


This    is   easily   illustrated   in    Figs.    3.15a   and 


hW(s) 


Hen(5) 


C 


W(s) 


Figure  3.15a 


Block  Diagram  of  the  Compensated  System 
(r=0)  . 
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Figure  3.15b.   Another  Configuration  of  Figure  3.15a 


In  order  to  obtain  the  elements  of  h,  solve  Eq. 

13 


3.66  or  equivalently ,  the  equation' 

l+G(s)H   (s)  =  [1  +  G(s)  G(-s) 
eq  R 


+ 


(3.70) 


where  the  [  ]   ( [  ]  )  indicates  the  poles  and  zeros  of 
1  +  G(s)  G(-s)  that  lie  in  the  left  (right)  half  s  plane. 

The  denominators  of  both  sides  of  Eq.  3.70  must 
be  equal  and  the  numerator  of  the  left-hand  side  can  be  ex- 
pressed in  terms  of  the  h.'s.  Therefore,  it  is  only  neces- 
sary to  find  the  L.H.P.  roots  of  [1  +  G(s)G(-s) ]  ,  either  by 

R 


13 


1+ 


G(s) 


vHT 


2  -  [i  +  G(s)G(-s) j+  fl  +  G(s)G(-s) j 


R  R 

From  Eqs.  3.70  and  3.71  it  is  possible  to  write 

[l+G(s)H   (s)][l+G(-s)H   (-S)]  =  [1+G(S)G?(~S)]  [1+G(s)G(-s)] 
eq  eq  k  k 

Assuming  that  G  (S)  has  poles  only  on  the  L.H.P.  it  can  be 

shown  that  the  poles  of  (l+G(s)H   (s) )  are  also  on  the  L.H.P 

eq 


and  therefore  l+G(s)H   (s)  =  [1  +  G(s)G(  s) ] 

eq  R 


+ 
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root-locus  or  factorization  techniques,  after  which  the  coef- 
ficients of  equal  powers  of  s  are  equated. 

The  results  obtained  by  Graham  in  two  given  prob- 
lems, by  this  method  and  Riccatti's  equation  method,  are  in 
very  good  agreement. 

D.   SHORT  COMMENTS  ON  OBSERVABILITY,  STABILITY,  CONTROLLABILITY, 
NONLINEAR  PROBLEMS,  STOCHASTIC  CONTROL  AND  PARAMETER 
IDENTIFICATION 

The  purpose  of  this  paragraph  is  to  give  some  references 
to  the  topics  enumerated  in  the  title. 

Because  Ref .  52  furnishes  all  the  necessary  Information 
about  the  references  prior  to  the  end  of  1969,  only  those 
published  after  this  date  will  be  considered. 

1.   About  Observability  and  Controllability 

Not  much  has  been  published  in  this  field,  and  so  far 
as  it  is  known,  all  the  reports  dealing  with  the  observability 
problem  are  already  mentioned  and  superficially  analyzed  in 
Ref.  52. 

Also  the  controllability  problems  have  been  in  a  peri- 
od of  abandon  until  recently  (19  70)  when  the  International 
Journal  of  Control  published  a  paper  by  Herget  [30],  written 
in  196  8.   The  author  considered  a  distributed  parameter  con- 
trol problem  of  a  linear  system  with  constrained  control, 
either  distributed  or  at  the  boundary. 

The  system  studied  is  represented  by  the  equation 

9Q(x,t) 

£t =  A(x)Q(x,t)  +  f(x)u(t)  ,   xefi,  teT  (3.71) 
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with  U(x)Q(x,t)  =  0   for  x  in  the  boundary,  and  Q(x,0)  =  0 
for  xcQ, . 

In  the  above  equations  A  (x)  represents  a  partial  dif- 
ferential operator  defined  relatively  to  Q,  U(x)  for  boundary 
x  is  an  operator  denoting  the  boundary  conditions  and  f(x)  is 
a  function  also  defined  in  Q. 

The  paper  shows  what  are  the  reachable  states  of  the 
above  type  of  system  for  two  classes  of  problems,  namely  with 
and  without  magnitude  constraint  in  the  norm  of  the  control 

(  |u||  2). 

2.   The  Stability  Problem 

In  the  field  of  stability,  some  recent  papers  seem  to 
be  of  great  interest  because  of  the  wide  variety  of  systems  to 
which  they  apply. 

The  following  references  are  mentioned  here: 
(i)   Kathri  [Ref.  35]  uses  multiple  Laplace  transform  tech- 
niques in  partial  differential  equations  which  can  be  written 
in  the  form 

I    ~i      K    J       ^k+j       M     ~m 

I    a.^v     j    £       9    v      £     3_v  =  (3>?2) 

i=0  13t1   k=l  j=l  k'^  8tk8x^   m=l  m  axm 
for  a  large  class  of  boundary  and  initial  conditions. 

The  system  was  reduced  to  a  transfer  function  expres- 
sion to  which  a  simple  stability  criteria  was  applied.   This 
report  is  quite  comprehensive  in  the  sense  that  the  author 
furnishes  a  step  by  step  explanation  of  all  the  details. 

Later,  the  same  author  published  a  new  paper  [Ref.  36] 
in  which  he  considered  the  same  type  of  systems  as  before  but 
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applied   to   sampled   data  problems   where    a  memory-less,    non- 
linear   feedback  was   present. 

(ii)       Ansari    [Ref.    2]    studied   the    sufficient    condition    for 
stability    of    a   system  described   by 

9QU,t)    +    u(t)3Q(x,t)    =    _   H(u)Tq(Xyt)     .      ,3-1,3,5 (3.73) 

for  0<x<L  ,  H>0  and  u>0  ,  where  the  control  u  =  f(T(L/t))  was 
such  that   q    and  Q (0 , t) f (Q (L,t) )  are  continuous  and  mono- 
tonically  decreasing  functions  of  u  and  T(L,t)  respectively. 
This  condition  can  be  verified  very  easily  and  although  ap- 
plicable only  to  a  small  class  of  p.d.e.'s,  englobes  some  im- 
portant diffusions  systems. 

(iii)   D'Souza  [Ref.  14]  treated  a  much  more  general  problem 
than  those  mentioned  previously.   He  considered  a  class  of 
causal  dynamic  systems  represented  by  nonlinear  vector-matrix 
p.d.e.'s.   For  a  good  understanding  of  the  developments  car- 
ried out  on  the  given  reference,  a  good  background  in  Func- 
tional Analysis  is  necessary,  with  special  emphasis  in 
Green's  function  type  of  problems.   Because  such  theory  is 
not  yet  within  the  range  of  the  knowledge  of  the  average  engi- 
neer, it  seems  to  be  useful  to  rewrite  the  paper  in  more 
practical  terms,  if  a  larger  audience  is  desired. 

(iv)   Kastenberg  [Ref.  37]  derived  the  stability  conditions 
for  nonlinear  feedback  control  systems  described  by  parabolic 
p.d.e.  of  the  following  type 

8Q(x,t)     M        X,    <    M       3Q(x,t) 

—it—  =  .1   aijM3a£  ax  +  I  bi(^^x—  +  c(x)+F(x,t,Q) 

(3.74) 
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where  F(x,t,Q)  represents  the  feedback  control  law,  which 
must  be  negative  definite.   The  initial  and  boundary  condi- 
tions are 

Q(x,0)   =   QQ(X)  for  x  inside  its  region  of  definition  and 
Q(x,t)   =   0      for  x  in  the  boundary  region  the  type  of 
system  studied  is  cf  wide  application  and  the  results  derived 
seem  to  be  easy  to  implement. 

3.  Nonlinear  Systems 

Although  published  in  1966 ,  the  paper  by  Uzgiris  and 
D'Souza,  "Optimal  Control  of  Distributed  Parameter  Systems 
with  Nonlinear  Boundary  Conditions"  [Ref.  6  7]  did  not  receive 
much  attention.   This  reason,  together  with  the  fact  that  the 
procedure  developed  is  relatively  simple  compared  with  the 
complexity  of  the  problems  it  can  solve,  motivated  the  present 
paragraph. 

The  authors  considered  the  case  of  a  one-dimensional 
linear  heat  equation  with  nonlinear  boundary  conditions  which 
were  made  to  include  the  bounded  control  function.   Then,  by 
discretizing  the  space  variable,  a  set  of  nonlinear  differen- 
tial equations  was  obtained  and  reduced  by  Laplace  transform 
and  convolution  techniques  to  a  nonlinear  vector  integral  equa- 
tion from  which  the  optimal  control  could  be  derived.   Follow- 
ing standard  methods,  the  final  problem  was  reduced  to  the 
solution  of  a  Volterra  type  equation  which  can  be  solved  by 
known  procedures . 

4 .  Stochastic  'Problems 

This  paragraph  considers  the  recent  papers  by 
Tzafestas  [Ref.  66]  and  Seinfeld,  Gavalas  and  Hwang  [Ref.  61] . 
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The  first  one  deals  with  the  identification  of  para- 
meters of  D.P.S.'s  by  reducing  the  system  to  a  lumped  para- 
meter system   for  which  there  exists  known  parameter 
identification  methods.   The  last  paper  contains  the  deriva- 
tion of  a  nonlinear  filter  for  systems  with  disturbances  in 
the  initial  and  boundary  conditions. 
5 .   Identification  of  Parameters 

The  present  chapter  is  concluded  by  mentioning  the 
modern  method  developed  by  Fairman  and  Shen  [Ref.  15]  in  the 
identification  of  D.P.S.   This  method,  named  by  the  authors 
as  the  "moment  functional  method"  is  applicable  to  one- 
dimensional  heat  or  wave  equations  and  in  the  case  of  the 
heat  equation  it  may  be  extended  to  the  problem  in  which  one 
coefficient  is  a  polynomial  function  of  time. 
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IV.   MODAL  CONTROL  THEORY 

The  theory  of  modal  control  was  first  developed  by 
Rosenbrock  [Ref.  53]  for  lumped  parameter  systems  and  the 
concept  generalized  by  Murray-Lasso  [Refs.  25  and  45],  and 
Gould  [Refs.  25  and  26]  in  order  to  include  D.P.S.'s. 

The  above  references  give  the  general  theory  of  what  is 
going  to  be  done  in  the  rest  of  this  thesis.   This  theory, 
although  applicable  in  many  situations,  suffers  from  several 
drawbacks.   The  major  ones  are  the  need  for  computing  the 
eigenvalues  and  eigenvectors  of  the  system's  operator,  calcu- 
lation of  the  inverse  of  singular  matrices,  need  for  r  con- 
trolling elements  if  r  modes  are  to  be  changed,  difficulty  of 
application  in  the  case  of  repeated  eigenvalues  and  also 
limitations  when  the  modes  are  complex  conjugates,  in  which 
case  only  the  real  parts  can  be  compensated. 

Some  of  the  inconvenient  parts  just  mentioned  were 
avoided  by  Simon  [Ref.  62]  and  Foster  [Ref.  18] ,  and  they 
furnish  very  interesting  research  material  in  the  optimal 
control  field. 

A.   GENERALITIES 

In  Fig.  1.1  is  shown  how  a  pointwise  input  produces  a 
distributed  output.   The  idea  of  the  modal  control  consists 
in  using  the  measurements  at  different  points,  uncoupling 
these  measurements  and  obtaining  the  coefficients  co.  of  the 
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14 
eigenf unction  expansion  of  the  output.     Then,  compensate 

these  eigenf unction ' s  coefficients,  compare  the  compensated 
values  with  the  reference  signal  and  build  with  this  informa- 
tion the  adequate  control  law.   This  is  done  in  a  synthesis 
procedure  conceptually  identical  to  the  analysis  procedure 
that  gave  the  coefficients  uj  .  .   The  global  picture  of  the 
control  system  for  a  heating  problem  is  indicated  in  Fig.  4.1 

and  it  shows  a  basic  assumption  that  was  made:   The  bandlim- 

15 
itedness    of  the  three  functions  considered  here,  namely 


14 

It  is  shown  in  Ref.  45  that  in  the  case  of  a  system 

described  by  an  operator  L  such  that  the  output  is 

y(x,t)  =  Lm(x,t)  and 

a)  L   exists, 

b)  The  eigenf unctions  of  L  are  products  of  functions 
of  time  and  distance, 

c)  L  is  completely  continuous  in  x,  that  is  to  say  it 
has  a  purely  discrete  spectrum  (eigenvalues  content)  and  the 
only  possible  cluster  point  is  the  origin,  then  it  can  be 
stated,  using  the  concept  of  generalized  Fourier  series,  that 


y(x,t)  =  I    o)i(t)ui(x) 


m 


(x,t)  =  I    Ui(t)ui(x) 


4-  V\ 

where  u. (x)  is  the  i    eigenfunction  of  L,  not  necessarily 
sine  or  cosine  functions,  and  forms  a  complete  set. 

Also  the  operator  L  is  such  that  its  adjoint,  L*,  has  a 
complete  set  of  eigenfunctions  v. (x) ,  which  is  orthogonal  to 
u.  (x)  . 

Using  the  orthogonality  properties  in  the  first  of  the 
above  equations  and  assuming  normalized  distance  and  eigen- 
functions it  follows  that 

1 
wi  (t)  =  /   y (x,t) v. (x)dx 
0 

15 

By  bandlimitedness  it  is  meant  that  the  functions  are 

considered  accurately  described  by  N  eigenfunctions,  instead 
of  the  infinite  number  that  exactly  characterizes  all  the 
partial  differential  equations. 
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reference  temperature,  output  temperature  and  control 
function. 


B.   LUMPED  PARAMETER  SYSTEMS 

The  basic  idea  in  the  modal  control  of  lumped  parameter 
systems  is  the  independent  shifting  of  the  lower  order  eigen- 
values such  that  the  speed  of  response  of  the  system  is 
increased. 

In  the  case  of  an  uncontrolled  linear  system  such  as 

x  =  A  x  (4.1) 

it  is  possible  to  obtain  x  as  a  linear  combination  of  terms 
containing  the  eigenvectors  and  eigenvalues  of  matrix  A.   For 
this  the  transformation 

y  =  H  x  (4.2) 

is  made,  from  which  it  follows  that 

(4.3) 

If  the  eigenvalues  of  A  are  distinct  it  is  possible  to 
obtain 


y    =   H   A  H    -1   y 


H    A   H         =    A   = 


\1     0 


'N 


and   therefore 


y .  =  a.  e 


X.t 

i 


(4.4) 


(4.5) 


where  a.  =  y .  (0)  . 

Finally,  from  Eq.  4.2  [Ref.  46,  p.  112] 

N        X.t 
x  =   y  a. u.  e 
-   nil  1~1 


(4.6) 
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Figure  4.1.   Theoretical  Implementation  of  the  Modal  Control 
of  a  D.P.S. 

Here  u.  is  the  i    eigenvector  of  A  and  \.    is  negative  if 
the  system  is  stable. 

The  purpose  of  the  modal  control  is  to  shift  the  eigen- 
values in  such  a  way  that  finally  the  following  relationship 
is  obtained 


H 


x  =   Y   a.  u.  ( 


(Xi+ki)t 


(4.7) 


with  k.  a  neqative  real  number  and  a.  the  distribution  of  the 
l      3  l 

initial  state  into' the  i    mode.   The  way  this  result  was  ob- 
tained is  explained  next. 
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1.   Ideal  Case 

The  system  is  described  by 

x  =  A  x  +  B  u  (4.8) 

and  the  output  is 

y  =  C  x  (4.9) 

In  this  ideal  case  the  dimensions  of  x  and  u  are  N 
and  the  dimensions  of  A  and  B  are  NxN. 

From  basic  Linear  Algebra  it  is  known  that  in  the 

case  of  square  matrices  with  simple  and  real  eigenvalues  the 

T 
eigenvectors  of  A  are  orthogonal  to  the  eigenvectors  of  A 

and  both  matrices  have  the  same  eigenvalues. 

From  the  above  statement  it  follows 


A  u.  =  X -u. 


A  v.  =  X  .  v. 


(4.10) 


and  representing  now  U  =  [u,  u  ,  .  .  .  uj  and  V  = 


Yi 


Y2 


Yn 


the  result  is 


A  U  =  U  A 

V  A  =  A  V. 


(4.11) 


When  the  different  u.'s  and  v.'s  are  normalized  it 

l        i 

can  be  written 

V  U  =  I  (4.12) 

and  from  this  equation  and  the  precedent  set  it  is  obtained 

A  =  U  A  V.  (4.13) 

Writing  now  the  state  equations  as  indicated  in 
Fig.  4.3 
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Figure    4.2.      Block   Representation   of    an   Uncontrolled   Linear 
System. 


A 


c 


K 


c 


Figure    4.3. 


Block  Representation  of  a  Controlled  Linear 
System. 
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x  =  (A  +  B  K  C)  x  (4.14) 

Finally,  because  of  Eq.  4.13,  a  convenient  choice  for 
B  and  C  would  be  making  them  equal  to  U  and  V.  which  qives 

x  =  U  (A  +  K)  V  x  (4.15) 

The  resultant  matrix  has  the  same  eigenvectors  as  A 
and  the  eigenvalues  have  been  increased  as  was  desired.   More- 
over, the  shifting  of  the  eigenvalues  is  realized  without 
interaction,  due  to  the  diagonal  form  of  A  and  K. 


16 


2.   The  Number  of  Manipulators  (r) 
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Is  Less  than  the 


Order  of  the  System  (N) 

Having  chosen  C  =  V  in  Eq.  4.14  the  same > equation  may 
be  written  as 

X  =  (UA+BK)  Vx  (4.16) 

Now,  in  order  to  have  B,  a  Nxr  matrix,  again  as  a 

z 

square    matrix,    make 


«-r-»- 


«-N-r* 


B    K    =     [U     i       0] 
~    ~  ~  r '       ~ 

Nxr  Nx(N-r) 


K 
~r 

i 

i 
i 

--L. 

0 

0 

1 
1 
I 
1 

5n- 

-r 

t 

r 
I 
t 
N-r 


(4.17) 


A   u   =    A    u 

Suppose  the  eigenvectors  are  different  and  call  them  w. 
Therefore, 

U(A  +  K)  V  w  =  (A  +  K)  w 

Aw   +   Kw   =    Aw   +   Kw 
But,    the   eigenvectors    of   A   are    unique    and   so  w   =    u    . 

17 

A  manipulator  is  the  device  to  which  the  control  law  is 
applied.   Its  output  is  the  physical  input  to  the  system  (e.g., 
heat  flow) . 
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Similarly, 


•r+   -<-N-r-> 


u  A  =  [u  :  a,  ] 

~  ..    L^r  « ~N-r 


~r 


i 
i 
i 
+. 

i 

i 

0   i 


V: 


r 

I 

t 

N-r 

4- 


(4.18) 


and  from  these  three  last  equations 


x  =  U 


A  +K  •  0 
r   r  i  ~ 


0 


£n- 


V  x 


(4.19) 


which   shows    how   it    is   possible    to   shift    arbitrarily   the    r 

lowest   eigenvalues. 

Reference    26,    pp.    246-24  8,    contains    the    theoretical 

development   for   the    situation   in   which   the    eigenvectors 

u.  , u      are   not   known    accurately. 

~i       ~r  2 


The  procedure  consists  in  making 

-1 


Br  =  D(VrD)  *  =  [bjbj br] 


(4.20) 


and  the  final  results  are  similar  to  the  ones  obtained  pre- 
viously, with  the  exception  that  disturbances  in  the  lower  r 
modes  cause  also  disturbances  in  the  higher  N-r  modes.   This 
is  not  always  a  serious  drawback  since  the  influence  of  the 
higher  modes  can  generally  be  neglected. 

3 .   The  Number  of  Sensors  (m)  is  Less  than  the  Order  of 
the  System  (N) 

The  derivation  for  this  case  is  also  in  Ref.  26.   It 

is  a  little  long  and  does  not  follow  exactly  the  same  lines 

as  the  distributed  parameter  case.   For  this  reason  this 

paragraph  will  only  be  concerned  with  the  basic  idea  of  the 

process . 
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Start  forming  an  Nxm  matrix 

?  =  [?1    '"    ?r    ~r+l  •*•  i?p  *?]_  •  '  *?ql       (4.21) 
where" p  is  the  total  number  of  modes  taken  in  account, 

q  =  m  -  p,  and  b  is  defined  as  in  Eq.  4.20. 

Then,  form  a  Nxm  matrix  E  such  that  e.  (i=l,...,p) 
approximates  as  closely  as  possible  of  u. (i=l , . . . ,p)  and 
of  b.  -  u.  for  the  last  q  elements  e.  of  E. 

Define  F,  nxm,  such  that 

FT  =  (ETW)~1ET  (4.22) 

This  matrix  may  be  written  as 


F  =  [f f   f  _,,..  .f   f .  .  .  .  f*] 

~i    ~r  ~r+l    ~p  -1    ~qJ 

=  [F  i  F    :  F    !  F*] 
q   r-q-  p-r'  q 


(4.23) 


and  from  it  obtain  C  defined  as 

*■  q  *    <-r-q->   f 

C   =  [F   +  F*  ';  F      ]n  (4.2  4) 

~q    q   r-q   + 

which  has  similar  properties  to  C  as  used  in  the  ideal  case. 
The  major  difference  is  that  the  rows  of  C  are  no  longer 
orthogonal  to  the  eigenvectors  u. (i=p+l, . . . ,N) ,  corresponding 
to  the  neglected  modes.   Therefore,  the  control  of  the  lower 
modes  interacts  with  the  higher  ones  not  allowing  a  perfect 
uncoupling.   Generally,  under  the  assumption  that  the  higher 
modes  are  sufficiently  high,  the  interaction  will  be  small. 
Figure  4.4  shows  the  block  diagram  illustrating  the 
present  situation.   Next  it  will  be  explained  how  the  fore- 
going concepts  can  be  extended  to  distributed  parameter 
systems . 
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Figure  4.4.   Block  Diagram  when  the  Number  of  Sensors  Is 
Less  than  the  Order  of  the  System 


C.   DISTRIBUTED  PARAMETER  SYSTEMS 

Before  starting  reading  this  section  one  should  become 
familiar  with  the  contents  of  Appendix  B. 

1.   Computation  of  the  Feedback  Control  Law 

Following  lines  similar  to  those  that  lead  to  Eq.  B.9 
it  results  that  in  the  case  of  m  and  y  being  both  functions 
of  x  and  t,  the  diagonalization  procedure  gives 


w.  (s)  =  X  .  (s)  y  .  (s) 


(4.25) 


where  oj  .(s)  and  y  .  (s)  are  the  Fourier  coefficients  of  y(x,s) 
and  m(x,s)  and  X  .  (s)  are  the  eigenvalues  of  the  operator. 

The  truncation  of  the  higher  order  modes  allows  to 
describe  the  control  system  in  matrix  form  as  shown  in  Fig. 
4.5. 

Calling  W  the  closed-loop  matrix  and  using  the 

properties  of  matrix  algebra  and  of  uncoupled  matrices  it 

can  be  written: 
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F 


c 


Figure  4.5.   Block  Diagram  of  the  Closed-Loop  System 


-1 


(4.26) 


W  =  (LF  +  I) 
From  this  it  follows 

(4.27) 
(4.28) 
and  finally 

W""1  =  F  +  L_1  (4.29) 

The  system  can  then  be  represented  as  in  Fig.  4.6  and 
for  each  loop  the  result  is 

(4-30) 


WL  -1  =  (LF  +  I)  1 
LW_1  =  LF  +  I 


Vs)  =  w-^iT  "  x-Ti7 


n 


n 


where  the  6  (s) ,  w  (s)  and  A  (s)  are,  respectively,  the 
n      n         n 

diagonal  elements  of  F,  W  and  L.   By  a  suitable  choice  of  the 
4>  's  it  is  therefore  possible  to  make  each  closed  loop  be- 
have as  fast  as  required. 

The  situation  illustrated  in  Fig.  4.6  although  desir- 
able is  not  possible  in  practice  because  of  the  existence  of 
sensors  restricted  to  certain  positions  and  of  manipulators 
which  cannot  give  exactly  the  desired  output  distribution.   It 
will  be  shown  next  that  in  a  physical  system  it  is  possible  to 
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Figure  4.6.   The  Uncoupled  System 

obtain,  the  coefficients  u).  from  the  output  measurements,  and 

a. ,  the  actual  input  to  each  manipulator,  from  the  knowledge 

of  the  distribution  produced  by  the  manipulators. 

2 .   Determination  of  the  Fourier  Coefficients  of  the 
Output  from  the  Output  Measurements 

Under  the  assumption  of  bandlimitedness 


N 


y(x,t)  =   Y  oo.  (t)  u.  (x) 
•  i    r      i 
i=l 


(4.31) 


Similarly,  the  measurements  at  x,,x2,...xN.  may  now  be  written 
as 


N 


y  (x  ,t)  =   J  oj.  (t)  u.  (x,) 
1       i±l  i      i   1 


(4.32) 


N 


y(x^t)  =   y  Ul(t)  ui(xN) 

i=l 
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or,    in   matrix    form 


Y  =  U  n 


(4.33) 


where 


Y   = 


yU-^t) 
y (x2,t) 


,     u 


ytx^t) 

It  follows  that 


u1(x1)u2  (x1)  .  .  .u.,  (x.^ 

U1(X2)  U2  (X2)  .  .  .VL-  (X2) 

U1(XN)U2(XN)"-UN(XN) 


,  n  = 


wx(t) 

032(t) 


o)N(t) 


(4.34) 


fi  =  U  XY 


(4.35) 
and  it  is  proved  in  Ref.  45,  Appendix  C,  that  there  is  always 

a  set  of  points  x,,x„,...x^  such  that  U  '  exists.   It  was 

verified  in  the  simulation  procedure  of  the  present  work  that 

even  in  the  case  of  the  sensors  being  relatively  close,  it  was 

possible  to  reconstruct  accurately  the  coefficients  go  .  .   The 

only  differences  are  in  having  higher  values  of  the  elements 

in  U    and  in  the  sensitivity  to  errors  in  the  measurements. 

There  are  several  criteria  described  in  Refs.  18  and  45  for 

the  minimization  of  the  sensitivity.   The  one  used  consists 

in  the  minimization  of  the  maximum  eigenvalue  of  the  product. 

(  U   ) X  U  (4.36) 

and  is  valid  whenever  the  number  of  sensors  (S)  is  greater 
or  equal  than  the  order  of  the  approximated  system  (N) . 
The  symbol  (I)  represents  a  pseudo-inverse  matrix  and  it  will 
be  defined  in  paradraph  d.   If  S=N,  the  pseudo-inverse 
coincides  with  the  inverse. 

When  S<N  the  sensors  must  be  positioned  according  to 


the  minimization  of  the  norm 


I  -  U   U  ( x . )   ,  whe  re  now  U 
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L 

is  function  only  of  one  element  of  the  set  {x.}»U   is  also 
a  different  matrix  and  it  can  be  obtained  as  explained  in 
Ref.  18,  pp.  60-6  3.   The  minimization  is  equivalent  to  ob- 
taining the  minimum  value  of  the  maximum  eigenvalue  of 
|I  -  U+  U  (x. ) I . 

3.   Determination  of  the  effect  of  the  manipulators 

Call  H  the  distribution  caused  by  the  i    manipulator 
and  decompose  it  in  an  eigenf unction  expansion. 
H-^x)  =  b11u1(x)+b12u2(x)  +  .  .  ,+b   ^(x) 

.*  (4.37) 

VX)  =  bNlUl(x)+bN2U2(x)  +  '-'+bNNUN(x) 
The  instantaneous  output  of  the  i    manipulator  (assuming 
instantaneous  response)  for  the  actual  input  a- (t)  is 
therefore 

a±(t)  Hi(x)  (4.38) 

and  given  that  the  desired  input  is 

N 
m(x,t)  =  I   u.(t)  u.  (x)  .  (4.39) 


Equating  b 

oth 

N 

I    «■ 
i-1  X 

(t)Hi(x)  = 

N 

j=i  D 

u. 

: 

,(x)  = 

N 

I  i 
i=l 

or 

T 
B  a 

= 

u 

and 

a  can 

be  obtained 

as 

a 

= 

(BT) 

-1 

N 


(4.40) 


(4.41) 


(4.42) 


T  -1  .  . 

The  problem  of  finding  (B  )    is  not  a  trivial 

one  and  for  many  types  of  manipulators  it  does  not  exist. 
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One  way  of  checking  for  the  existence  of  the  inverse  is  to 
compute  the  determinant  of  B;  if  its  value  is  very  small  it 
is  better  to  try  to  avoid  it  because  a  tremendously  high  con- 
trol gain  will  be  required. 

By  far  the  best  situation  is  when  the  distributions  of 
the  manipulators  coincide  with  the  eigenf unctions ,  in  which 
case  B  will  be  the  identity  matrix.   This  situation  is  not 
very  realistic;  however,  it  is  frequently  possible  to  move  the 
position  of  the  manipulators.   In  this  case,  the  best  obtain- 
able positioning  (n • )  turns  out  to  be  the  one  which  maximizes 
the  projection  of  H.fx-n.)  on  u. (x) . 

The  adequate  choice  of  the  manipulators  is  also  a 
question  of  good  sense,  and  before  starting  to  move  them  and 
trying  to  maximize  the  above  mentioned  projection,  it  is  high- 
ly desirable  to  look  physically  at  the  distribution  and  see  if 
there  is  any  possibility  of  obtaining  it  with  the  available 
manipulators.   As  an  example,  it  may  be  said  that  given  a  sys- 
tem described  by  Eq.  2.7  it  would  be  absurd  to  try  to  get  the 
heat  configuration  shown  in  Fig.  4.7  with  the  heat  source 
positioned  as  shown. 

It  is  now  possible  to  represent  the  control  of  the 
given  system  in  an  easily  implementable  way,  as  illustrated 
in  Fig.  4.8. 

It  should  be  mentioned  that  the  coefficients  t .  are 

not  the  Fourier  coefficients  of  the  desired  output  distribu- 

1 
tion  (t  .  (t)  =  /  i (x,t) v .  (x) dx) )  ,  but  these  values  multiplied 
3  0         : 
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1   X 


Figure  4.7.   Example  of  an  Undesired  Control 

1+4)  A  G 

by  — ,  J nJ — J ,  the  inverse  of  the  closed-loop  gain; 

D  D 

G.  represents  the  gain  of  the  j    manipulator.   This  statement 
is  also  not  always  exact  as  will  be  seen  later  in  the  analysis 
of  the  computer  results.   Meanwhile  it  may  be  added  that  in 
many  cases  a  minimum  square  error  deviation  may  not  be  de- 
sired, but  another  kind  of  criteria;  in  such  conditions  the 
t's  must  be  chosen  using  an  optimizing  program,  as  for  example 
the  gradient  search.   When  the  number  of  filters  is  too  high 
this  search  probably  will  require  a  great  deal  of  computer 
time  and  it  may  become  difficult  to  determine  if  the  obtained 
optimum  is  local  or  absolute.   This  procedure  was  implemented 
successfully  for  three  reference  coefficients  and  one  hundred 


18 

The  product  of  the  coefficients  x .  by  the  inverse  of 

thesis 


loop  gain  will  be  mentioned  throughout  the  rest  of  this 
as  the  "compensated  Fourier  coefficients." 
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solutions  of  the  system  only  took  two  minutes  of  computer 
time . 

4 .   Practical  Limitations 

As  in  the  lumped  parameter  case  the  most  common  situa- 
tion takes  place  when  the  order  of  the  approximated  system 
differs  from  the  number  of  sensors  or  manipulators.   In  this 
case  the  problem  of  inverting  non-square  matrices  must  be  con- 
sidered and  it  can  be  solved  using  the  concept  of  the  pseudo- 
inverse  (  +  )  . 

The  rules  for  calculation  of  the  pseudo  inverse  of  real 
matrices  are  simply  synthesized  as  follows.   For  further  in- 
formation see  Refs.  42  and  73. 

4- 

a.  If   A   is    in    diagonal    form,    A      is    also   diagonal,    with    the 

non-zero  elements    the . reciprocal    of   those    of   A   and  with 
zeros   where    A  has    zeros. 

m  i 

b.  Knowing  that  A  A  is  a  square  matrix  obtain  A   as 

A+  =  (ATA)"1AT  =  AT(AAT)_1  (4.43) 

The  computer  implementation  of  these  concepts  is  shown  in 
the  next  chapter.   Depending  on  the  numerical  values  of  the 
elements  in  the  matrices  sometimes  one  procedure  is  better 
than  the  other  one.  * 
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V.   DETAILED  SOLUTION  OF  A  PROBLEM  BY  MODAL  CONTROL  TECHNIQUES 

A.   PROBLEM  DESCRIPTION 

It  is  desired  to  heat  a  given  rod  at  three  different 
points  such  that  a  certain  temperature  distribution  is 
achieved. 

Geometrically/  Fig.  5.1  shows  the  system  configuration. 
Looking  at  one  differential  volume  element  it  is  possible  to 
derive  the  equations  describing  the  process. 


Heat  Input  (q0^,t)) 


x  =  0 


x=1 


Figure  5.1.   The  Heating  of  a  Rod 

qa(x„t) 
q2(x„t) 
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dS    - 

dV    - 


*!  dx, 


temperature 
heat      flow 
cross    section 
lateral   element 

of   area 
volume    element 

of    length    dx 


Figure  5.2.   Heat  Flux  in  a  Volume  Element 


Consider  the  principle  of  conservation  of  energy;  from 
Fig.  5.2  the  following  equation  can  be  written: 

9E 


(qrq2)A  +  q3dS  =  ^dV 


(5.1) 


10  3 


However, 


3q, (x    ) 
q2(x1)    =    qi(x1+dx1)     =   q1(x±)    +    — ±-dXl  (5.2) 

and  this    gives 

3q2 

qrq2  "  -  ■^dxi  <5-3> 

Taking  E=cpy.,  where  c  is  the  specific  heat  at  constant 
volume  and  p  the  mass  density  it  follows 

9qx  3yx 

"  8xY  Adxl  +  q3S  dxl   =   CP^A  dxl  (5'4) 

or 

a      8yl    2  8yl    ? 

-  g~-  (-k  3^-)ttR  dxx  +  q3  2-rTRdx1  =  cp~  ttR  dx±  (5.5) 

where  k  is  the  thermal  conductivity  (Btu/sec-ft-°F) .   Elimi- 
nating the  common  terms  and  taking  a  =  k/pc  (a  is  the  thermal 
dif fusivity) : 

^2Y±        2q3    3y1 
a 7T   +    — =■  =  7TT—  (5.6) 

8x2    pcR    dt]_ 

This  equation  can  be  normalized  and  simulated,  taking  into 
account  the  scaling  factors,  as 

M1  +  q  =  !¥■  ,       (0<x<l,  t>0)  (5.7) 

8x2        dt 

For  simplicity  the  following  initial  and  boundary  condi- 
tions are  considered: 

y(x,0)  =  0  ,   0<x<l  (5.8) 

y(0,t)  =  0 
a  }   0<t<«>  (5.9) 
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Assume  the  availability  of  two  different  sets  of  three 
manipulators,  one  giving  a  heat  distribution  coincident  with 
the  eigenf unctions ,  as  shown  in  Fig.  5.3,  and  the  other  giving 
the  heat  distribution  illustrated  in  Fig.  5.4. 


fl  X 


Figure  5.3.   Heat  Distribution  No.  1 


Heat 


1  x 

Figure  5.4.   Heat  Distribution  No.  2 

Throughout  the  simulation  runs  only  the  two  desired  tem- 
perature distributions  of  Figs.  5.5  and  5.6  will  be  considered, 
the  first  one  of  which  has  the  shape  of  the  first  eigen- 
f unction  and  being  for  this  reason  easily  attainable. 
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Temp 


1  x 


Figure  5.5.   Temperature  Distribution  No.  1 


Figure  5.6.   Temperature  Distribution  No.  2 

B.   DERIVATION  OF  THE  EIGENFUNCTIONS  AND  EIGENVALUES 

Consider  again  Fig.  B.l.   By  Laplace  transforming  the 
distance  dependent  operator  with  respect  to  t  it  follows  that 
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m(x,t) 


Figure  B.l.   Representation  of  the  System  in  Operator  Form 


M[-]  =  L_1[.]  =  < 


[—^2  +  sH-> 

dx         ,  .    . 

T(0,s)  =  0;  dT;X'SJ  1=  0 

dx    '   . 
x=l 


To  find  the  eigenfunctions  make 

M  u  =  X  u 


or 


(         2 
d_u 

dx 


+  (s-X)u=0 


V 


u(0,s)=0,  du(^l|    = 

x=l 


whose  solution  is 


(5.10) 


(5.11) 


(5.12) 


j/X^ 


-j/X  -sx 


From  the  boundary  condition  u(0,s)=0: 


(5.13) 


u(x,s)  =  C  sin/X  -  s  x 
and,  from  the  other  boundary  condition, 

A    -s  =  (2n+l)5- 


n 


or 


X   =  s  +  (2n+l)2^- 
n  4 


(5.14) 


(5.15) 


(5.16) 


In  order  to  normalize  the  eigenfunctions  take  <u  ,u  >  =  1  or 

3  n '  n 


/  u  (x,s)  dx  =  1 .   The  following  expression  results: 
0   n 
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C2  /  sir^/X^s  x  dx  =  1  ,  (5.17) 

from  which  C  =  /2 
n 

The  next  step  is  the  determination  of  the  adjoint  operator, 
M  ,  by  using  the  common  technique  indicated  in  Ref.  45,  pp.  269' 
271  and  by  Crandall  [Ref.  12,  p.  211]. 

By  definition  of  adjoint  operator  write 

<M  u,  v>  =  <u,  M  v>  (5.18) 

Integrating  by  parts  the  left-hand  side: 

r    r      «  u    .     ,„   ,  %     !  -,  du      ,  ,      du   dv    , 

/     [ j  +    (s-^)ul*v   dx   =    "      d3Tv]       +    /      dx~'dx~  dx   + 

0  dx  ax        0         0 

(5.19) 

+  /    (B-X)-vdx  :=    [-   g-v]    +u      &       -    /1u.^4+    /1(s-X).vdx 
0  ax        0        L  ax   0        0        dx2        0 

From  the   equations    just   above    the    definition    of    adjoint 

operator   is    verified   taking 

+  d2 

M      = =*-  +    (s-X)  (5.20) 

dx^ 
and 

v(0)    =    — |         =0  (5.21) 

dX    x=0 

+  + 

Because  M  =  M  and  the  boundary  conditions  for  M   are  the 

same  as  those  for  M,  the  former  is  said  to  be  self-adjoint. 

This  fact  permits  writing: 

1 
M[«]  =  I   X  (s)  u  (x)/  v  (O         dc  (5.22) 

111.  II       <-.     II 

0 

and,  according  to  an  important  property  of  the  functions  of 
an  operator  (see  also  Eq.  Bll) . 
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1 

f(M[-])  =  I    f(\  (s))u  (x)/  v  (C)[*]d  c  (5.23) 

n=l    n     n    0   n 

which  means  that  if  f(M[«])  =  M   [•]  =  L[«] 
it  follows  that 

Ll-]  -  I,    FWUn(x)/  vn(C)  [-]  d  ?     <5-24) 

n=l  n        0 

Next,  compute  the  first  six  eigenvalues  from  Eq.  5.17: 

i       A  .        /T7=^" 
i  i 


0  s+2.4674  1.57080 

1  s+22.2066  4.71239 

2  s+61.6850  7.85398  (5.25) 

3  s+120.9027  10.99558 

4  s+199.8595  14.13717 

5  s+298.5555  17.27876 

The  corresponding  first  three  eigenvalues  of  the  operator 
L  are 

Xo(s)   =   s+2.4674 

X1(S)   =   s -t-22. 2066  (5'26) 


X0(s)   = 


1 


s+61.6850 

C.   COMPUTATION  OF  THE  FEEDBACK  CONTROL  LAW 

Suppose  it  is  desired  to  make  the  first  three  eigenvalues 

coincide  with  the  fourth  one.   Using  Eq.  4.30  it  follows: 

$  =  s  +  120.9027  -  s  -  2.4674  =  118.4353 
o 

4>1   =  s  +  120.9027  -  s  -  22.2066  =  98.6961      (5.27 
<J>   =  s  +  120.9027  -  s  -  61.6850  =  59.2177 
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D.       COMPUTATION    OF   THE    SENSORS'    POSITIONS    AND    OF    U_1 

Consider   now   the    cases    of   having    three    and   six  sensors. 
The    second   of  Eqs.    4.35    is    repeated 

u1(x1)    u2(x±)    - u^) 

ul(x2]    u2(x2}     UN(X2) 


U   = 


U;L(xs)    u2(xs) 


UN(XS) 


(4.49) 


where  x,,x_,...,xs  are  the  sensors'  locations. 

Two  computer  programs  were  written  for  the  conditions 
S>N.   The  first  one  divides  the  interval  x(0,l)  in  ten  parts 
and  tries  all  the  possible  combinations  of  x~  ,x~ ,x~   within 
this  interval,  such  that  0<x,  <x«<x^<_l .   If  the  temperature 
at  x=0  was  not  known,  the  possibility  of  a  sensor  there  should 
also  be  considered.   Once . the  optimal  rough  positions  are  com- 
puted, this  information  is  introduced  in  a  gradient  subroutine 
which  searches  for  the  rigorous  optimal  positions.   Both  pro- 
grams make  use  of  the  same  function  subprogram  (EVMAX) ,  which 
in  turn  uses  two  IBM  library  subroutines  (MINV  and  MPRD, 
respectively  for  the  inverse  and  product  of  matrices)  and  one 
N.P.G.S.'s  subroutine  (JACVAT)  for  the  computation  of  the 
eigenvalues  of  a  real -symmetric  matrix. 

The  optimal  positions  and  corresponding  U  and  U   matrices 
derived  for  the  cases  of  three  and  six  sensors  are: 


CASE  I- 


x(l)  to  x(3)   =   .28570E00   .57144E00    .85715E00 

(5.28) 
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Ull  t0  U13   =  -61357E00    .13787E01    .11058E01 

U21  t0  U23   =  •11057E01    .61353E00   -.13788E01     (5.29) 

U31  t0  U33   =  •13788E01   -.11057E01    .61367E00 

Ull  t0  U13   =  -17532E00    .31593E00    .39391E00 

U21  t0  U23   =  -39393E00    .17533E00   -.31591E00     (5.30) 

U3J  to  U33   =  .31589E00   -.39391E00    .17532E00 

CASE  II- 

x(l)  to  x(3)  =   .15406E00    .30828E00    .46234E00 

x(4)  to  x(6)  =   .61621E00    .76996E00    .92325E00 


(5.31) 


Ull  t0  U13  =  -338904E00  .938860E00  .132315E01 

U21  t0  U23  =  -658373E00  .140436E01  .932899E00 

U31  t0  U33  =  -939128E00  .116083E01  -.665084E00  (5.32) 

U„,  to  U„0  =  .116491E01  .333099E00  -.140276E01 
41      43 

Uc,  to  Uco  =  .132288E01  -.661515E00  -.330576E00 

U.,  to  U,^  =  .140394E01  -.132272E01  .116496E01 
61      6  3 

U*   to  U*  =  .05225   .10148   .14464  .17915   .20310   .21531 


11     "16 

4-       4- 

J21  t0  U26 

4-         I 

'31  t0  U36 


U21  t0  U2f 

U*   to  U*    =   .20383   .14385  -.10216   .21588  -.05109   .17894 


(5.33) 

E.   COMPUTATION  OF  THE  REFERENCE  COEFFICIENTS 

Using  the  definition  of  x .  write 
y  3 

1 

t  .  (t)  =  /  i(x,t)v.  (x)dx  (5.34) 

3  0         3 

where  i(x,t)  is  the  desired  temperature  distribution.   This 
gives,  for  the  two  desired  temperature  distributions  shown  in 
Figs.  5.5  and  5.6/  respectively 
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TAU  (1)  =  v-j     =      0.707107 

TAU  (2)  =  0.0  (5.35) 

TAU  (3)  =  0.0 
and 

TAU  (1)  =  0.960317 

TAU  (2)  =  0.685931  (5.36) 

TAU  (3)  =  -0.32011 

F.       DERIVATION    OF    B    FROM   THE    HEAT    DISTRIBUTION    OF    THE 
MANIPULATORS       ~ 

For   the   heat   distribution    coincident  with   the   eigenfunc- 
tions,    B    is    the    identity  matrix. 

In   the    case    of   the    heat    distribution  No.    2,    obtain    from 

Eq.    4.38    (using   the    orthonormality   properties)    and   Fig.    5.4 

the    coefficients    of   B    as    follows: 

0.5 
b  =      /2    /        sin    2ttx    sin^dx      =      .169765 

11  0 

0.5  3tt 

b,  9      =      /2    j         sin    2ttx   sin-^-dx      =       .363783 
LZ  0 

0.5  - 

b  =      /2~  /        sin    2ttx   sin^dx      =       .353678 

1J  0 

1 
b  =      /2~  /    sin    ttx   sin^dx      =      .600211  (5.37) 

0 

1  . 

b  =      fZ  j    sin    ttx    sin^dx      =       .360126 

11  0 

1  , 

b  =      /2~   /    sin    ttx   sin^dx      =      -.085744 

16  0 

1 


b-,       =      -/I   /    cos    ttx   sin    ^dx      =       .424413 


1 
0 


'32      =      -/Z   j    cos    ttx   sin    ^dx      =      -.381972 


1  5TT 

b_~      =      -/2~  /    cos    ttx   sin   ^dx      =       .151576 
3  0 
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T  -1 
After  this  Q  =  (B  )    must  be  computed.   The  result  came 


out  to  be  the  following  one 


0(1,1)  to  0(1,3) 
Q(2,l)  to  Q(2,3) 
Q(3,l)  to  0(3,3) 


-.12283 

1.07012 

.89194 


.71648 

.69963 

-1.27602 


2.14944 

-1.23328 

.88434 


(5.38) 


G.   NUMERICAL  SOLUTION  OF  THE  PARABOLIC  DIFFERENTIAL  EQUATION 
THAT  DESCRIBES  THE  SYSTEM 

The  system  was  simulated  using  the  Crank -Ni col son  method. 

The  following  system  of  equations  resulted,  where  Jl  stands 

for  the  instant  of  time  J+l.   The  normalized  distance  interval 

was  divided  in  20  parts  and  the  time  interval  was  made  equal 

to  1/400. 


4-1 
-1  4-1 
-1  4-1 


•1  4-1 
-2  4 


u 


2,  jl 


u 


u 


3,jl 


4,jl 


u 


If  jl 

0 

0 


u 


20,  jl 


u 


21, jl 


u 


Ifj 


u 


u 


2,j 


3,j 


h 


u 


19, j 


u 


3,j 


4,j 


u 


5,j 


2u 


20, j 


l2,j 


3,j 


[4,j 


+26t 


u 


21, j 
0 


l20,j 


L21,j 


(5.39) 


This  system  is  solved  in  double  precision  by  IBM  sub- 
routine "DGELB."   In  order  to  take  in  account  the  fact  that 
the  sensors  are  not  generally  in  the  exact  subdivisions  of 
the  distance  domain,  also  IBM  subroutine  "ALI"  was  used,  which 
interpolates  for  the  sensors1  positions. 
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In  Eq.  5.39  the  letter  u  stands  for  temperature  (contrar- 
ily  to  the  rest  of  the  chapter  where  the  temperature  is 
represented  by  y)  in  order  to  be  consistent  with  the  computer 
program. 

H.   SIMULATION  RESULTS  BASED  IN  KNOWN  THEORETICAL  FACTS 

Ten  different  simulations  were  performed.   These  simula- 
tions correspond  to  changes  in  the  following  parameters: 
feedback  gains,  manipulators  output  shape  and  gain,  reference 
signal  and  number  of  sensors.   The  runs  involving  the  use  of 
strong  gain  saturation  and  a  control  different  from  the  eigen- 
f unction  control  gave  origin  to  the  development  of  new  prac- 
tical procedures.   Due  to  the  relevancy  of  such  a  fact,  the 
next  chapter  will  be  entirely  devoted  to  it. 

This  thesis  shows,  so  far  as  it  is  known,  the  first  closed- 
loop  simulation  of  the  modal  control  of  a  distributed  parameter 
system  by  non-optimal  control  techniques.   In  order  to  illus- 
trate the  characteristics  of  this  model  the  present  section  is 
dedicated  to  emphasize  the  results  of  changes  in  the  three 
most  representative  parameters,  namely  the  feedback  gains, 
the  direct  gains  and  the  reference  coefficients. 

With  exception  of  run  No.  7  all  the  other  ones  were  made 
using  the  Crank-Nicolson  method,  as  described  in  Chapter  II. C. 
1.   This  is  a  brute  force  method  that  has  the  advantage  that 
it  is  possible  to  make  it  quite  accurate  only  by  reducing  the 
time  and  distance  intervals,  but,  in  turn,  it  requires  a 
sophisticated  computer  program  and,  consequently,  a  long  com- 
puter time.   The  program  written  makes  use  of  I.B.M.  subroutine 
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"DGELB"  for  the  double  precision  solution  of  a  system  of 
linear  equations,  and  "ALI"  which  takes  the  temperature  at 
the  descritized  points  in  the  rod  and  interpolates  these  val- 
ues in  order  to  get  the  temperatures  at  the  sensors1  positions. 
At  the  end  of  this  thesis  is  shown  the  computer  program  used 
in  run  No.  10,  which  serves  as  reference  program  for  the  other 
runs . 

Run  No.  7  was  done  using  the  Bubnov-Galerkin  transformation 
and  this  permitted  writing  a  fairly  simple  program  using  a 
prediction-correction  numerical  solution  of  ordinary  differen- 
tial equations  [Milne,  Ref.  44].   This  technique  can  very 
easily  be  applied  to  a  large  variety  of  systems,  as  described 
by  Foster  [18] .   However,  when  dealing  with  the  modal  control 
approach,  the  matrix  B  of  the  Bubnov-Galerkin  method  cancels 
the  matrix  Q=B    in  the  feedback  loop.   This  makes  the  trans- 
formation only  applicable  in  the  case  of  having  the  manipulator 
distribution  coincident  with  the  eigenfunctions ,  for  which  B  is 
the  identity  matrix.   For  this  reason  Foster  uses  the  method 
only  in  an  optimal  control  way,  applied  to  the  linear  regu- 
lator problem,  such  that  the  mentioned  transformation  does 
not  take  place. 

1.   Run  No.  1;   Feedback  Gain  Equal  to  Unity 

The  differences  from  the  reference  program  are: 

a.  Dimension  -  U(21,200) 

b.  Data  -  NT  =  200 

c.  TAU(l)  =  .707107 

d.  PHI(l)  =  PHI(2)  =  PHI  (3)  =  1.0 
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e.  DO  38  10=1,  Nl 

EVl(ID)  =  1.57080*H*(ID-1) 
EV2(ID)  =  4.71239*H*(ID-1) 
EV3(ID)  =  7.85398*H*(ID-1) 
Hl(ID)  =  1.414214*SIN(EV1(ID) ) 
H2(ID)  =  1.414214*SIN(EV2(ID) ) 
38H3(ID)  =  1.414214*SIN(EV3(ID) ) 

f.  First  IQ  is  4,  second  IQ  is  40  and  the  last  is  200 

g.  Q  is  the  identity  matrix 

h.   The  loop  "DO  41  ..."  was  not  included  and  this 
allowed  negative  control. 

In  this  run  it  was  intended  to  obtain  the  temperature 
distribution  No.  1  (Fig.  5.5)  using  manipulator  control  No.  1 
(Fig.  5.3);  the  reference,  temperature  coefficients  were  not 
divided  by  the  respective  closed-loop  gains.   As  it  can  be 
observed  the  output  never  approached  the  desired  distribution, 
converging  very  slowly  toward  an  amplitude  much  lower  than  the 
desired  one;  this  illustrates  the  need  for  using  the  compen- 
sated Fourier  coefficients  as  described  in  Chapter  IV. 3. 

2 .   Run  No.  2:   Increasing  the  Feedback  Gain  of  the  First 
Ei gen function 

Same  statements  as  in  previous  run  except  for  PHI(l)  = 

5.0.   The  convergence  although  faster  than  in  the  preceding 

case  was  still  too  slow.   Notice  that  the  closed-loop  gain 

corresponding  to  the  first  eigenfunction  decreased.   This  run 

illustrates  how  the  total  gain  of  the  system  decreases  when 

the  feedback  is  increased. 


3  16 


a)      T  =   0.01 


1         X 


b)       T    =    0.1 


c)       T    =    0.5 


1       X 


Figure    5.7.       Run   No.     1 
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a)       T   =    0.01 


1         X 


b)       T    =    0.1 


1      X 


c)       T    =    0.5 


1      x 


Figure    5.8.       Run  No.    2 
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3.  Run  No.  3:   Increase  of  the  Manipulators'  Gains 

The  differences  relative  to  case  No.  1  are  as  follows: 

c.  TAU(l)  =  3.62272 

d.  PHI(l)  =  5.0 

f.   First  IQ  is  1,  second  IQ  is  2  and  the  third  one  is 
40. 

Also  in  statement  4  3  the  total  value  was  multiplied  by 
20,  which  is  equivalent  to  an  increase  of  the  manipulator's 
gains  by  the  same  factor.   The  response  was  very  fast  and  at 
T=.02  sees  the  system  was  practically  stabilized  at  the  de- 
sired temperature  distribution. 

The  indicated  value  of  TAU(l)  was  found  after  having 
divided  the  original  TAU(l)  by  the  respective  closed-loop 
gain.   The  result  of  this  division  is  3. 62277,  practically 
equal  to  the  3.62272  that  gave  the  exact  desired  heat  distri- 
bution.  The  reason  for  this  small  difference  is  due  to  the 
fact  that  the  matrices  P  and  Q  are  approximate  transforma- 
tions (in  the  specific  case  of  eigenfunction  control  consider 
only  the  effects  of  P)  and  maybe  due  to  round-off  errors.   In 
general,  if  N  is  high  the  empirical  and  theoretical  values  must 
be  very  close  to  each  other. 

4 .  Run  No.  4:   Making  the  First  Three  Eigenvalues  to 
Coincide  with  the  Fourth  One 

The  statements  below  are  the  only  changes  in  the  last 


run. 


PHI(l)  =  118.5353 

PHI (2)  =  98.6961 

PHI (3)  =  59.2177 

TAU(l)  =  85.502  (derived  value  was  85.503) 
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a)       T   =    0.01 


1        X 


b)       T  =    0.1 


1         X 


T    * 
1 


WANT -HAVE 


c)       T   =    0.5 


0 


1    x 


Figure    5.9.       Run   No.    3 
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a)       T   =    0.0025 


1       X 


b)       T   =    0.005 


1       X 


c)       T    =    0.1 


1       X 


Figure    5.10.       Run   No.    4 
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Again  statement  4  3  has  the  normal  gain.   The  results  were 
practically  the  same  as  in  the  case  just  considered,  with  the 
advantage  that  much  less  heat  power  was  required.   When  trying 
to  increase  the  gain  of  the  manipulators  by  a  factor  of  20  the 
system  became  unstable. 

Some  other  modifications  were  done,  as  for  example 
PHI (2) =PHI (3) =0  but  no  reasonable  change  was  observed,  which 
shows  the  dominant  influence  of  the  first  eigenfunction  in 
this  specific  problem;  it  is  not  so  in  many  other  cases.   An- 
other change  was  in  the  position  of  the  sensors;  for  x(l)= 
.3,x(2)=.4  and  x(3)=.5  a  different  corresponding  matrix  P 
with  much  greater  values  was  obtained.   As  expected,  the 
transformation  carried  on  by  P  gave  the  necessary  eigenfunc- 
tion information  and  the  results  are  practically  the  same  as 
when  the  sensors  were  in  the  optimal  positions  (they  differ 
only  in  the  fourth  decimal  place) .   The  important  difference 
between  these  cases  resides  in  the  fact  that  the  optimal  posi- 
tions are  the  ones  for  which  the  errors  in  the  measurements 
are  minimized. 

5 .   Run  No.  5:   A  different  Output  Distribution  and 
Nonlinear  Control 

In  order  to  have  temperature  distribution  No.  2  (Fig. 

5.6)  the  last  run  was  modified  as  follows: 

90  DO  91  IH=1,N1 

X1(IH)=(IH-1) *H 

D=6.283185*X1(IH) 

Y2(IH)=1.0-COS(D) 

91   CONTINUE 
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a)       T   =    0.0025 


1       X 


b)       T   =    0.005 


1     X 


C)       T    =    0.1 


1        X 


Figure    5.11.       Run   No.    5 
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A   saturated    amplifier  with    -20<_a<_20   was    used.       It   re- 
sulted  that    although  this    control    is   quite    saturated,    a   good 
output"    (at    steady-state    it  equals    the    one    obtained,    with    linear 
control)    was    obtained  with   the    analytically    derived   x's, 
namely: 

TAU(l)    =    116.10 

TAU(2)    =       82.93 

TAU(3)  =  -38.70 
Statement  40  was  substituted  by  the  sequence: 

40  DO  41  1=1, MP 

IF  (ALPHA(I) .GT.20.0)  ALPHA (I) =20 . 0 

41  IF  (ALPHA(I) .LT.-20.0) ALPHA (I) =-20 . 0 
6 .   Run  No.  6:   Larger  Number  of  Sensors 

The  number  of  sensors  was  increased  from  three  to  six. 
This  required,  relative  to  the  preceding  case,  the  following 
changes  in  the  computer  program: 

DIMENSION  -  x(6)  ,P  (3,6)  ,TEMP(6) 

DATA       -  KS=6 

x(l)  to  x(6)  as  in  Eq.  5.31  and  P  as  in  Eq .  5.33 

TAU(1)=116.10 

TAU(2)=  82.9  3 

TAU(3)=-38.70 

The  loop  40  DO  41  etc.,  was  removed  and  again 

statement  40  became  the  same  as  in  the  reference  program. 

This  run  was  done  in  order  to  illustrate  the  technique 
of  using  the  pseudo-inverse  matrix. 
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WANT 


a)    T   =    0.0025 


1        X 


WANT 


b)       T   =    0.005 


1       X 


c)       T   =   0.1 


Figure    5.12.       Run   No.    6 
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7 .   Run  No.  7:   The  Bubnov-Galerkin  Transformation 

The  computer  program  for  this  case  is  also  shown  in 
the  end  of  the  thesis.   The  loop  elements  are  exactly  as  in 
paragraph  five  but  the  nonlinearity  was  removed.   Practical- 
ly no  difference  is  noticed  in  steady  state,  which  is  far 
better  than  expected  and  shows  that  in  the  present  example 
the  higher  order  modes  of  the  system  can  be  neglected. 
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a)       T   =    0.0025 


b)       T   =    0.005 


c)       T   =    0.1 


1      X 


Figure    5.13.       Run   No.    7 
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'-  VI  .   SIMULATION  RESULTS  BASED  ON  A  NEW  TECHNIQUE 

It  is  a  mathematical  fact  that  the  Fourier  coefficients 
of  an  N  terms  expansion  minimize  the  square  error  of  the  ap- 
proximation.  When  the  matrices  P  and  Q  are  high  order 
matrices  (theoretically  infinite  order)  the  use  of  the  compen- 
sated Fourier  coefficients  for  input  reference  gives  exactly 
the  desired  output.   This  may  be  the  reason  why  Gould  states 
[Ref.  26,  p.  242]  that  typically  the  number  of  sensors  is  of 
the  order  of  fifty  and  the  number  of  manipulators  of  the  order 
of  ten.   When  the  analysis  described  by  P  is  already  relatively 
accurate  (this  happens  in  the  example  given  in  this  thesis  even 
with  only  a  third  order  matrix)  and  the  synthesis  done  by  Q  is 
also  precise  (is  exact  for  eigenfunction  control)  the  compen- 
sated Fourier  coefficients  also  give  origin  to  an  output 
distribution  close  to  the  desired  one.   For  a  better  under- 
standing of  what  follows  the  above  examples  are  nemed  as  ideal 
cases.   However,  when  the  accuracy  of  the  transformations 
defined  by  P  and  Q  is  inadequate  it  is  possible  to  improve  it 
by  changing  the  reference  coefficients  under  the  control  of  a 
gradient  search  until  an  output  distribution  that  better 
approaches  the  desired  one  is  achieved.   In  order  to  be  able 
to  use  the  gradient  search,  similarly  to  what  was  done  for 
the  computation  of  the  optimal  sensors'  positions,  a  certain 
cost-function  (SQUERR)  was  defined. 
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If   it   is    intended   to   minimize    the    square    error  between 
the    desired    and    obtained   distributions,    one    convenient    defini 
tion.of    SQUERR   is,    according   to    Fig.    5.14: 

M 

SQUERR      =  Y     (yl.     -    y2.)Z 

The  computer  program  used  to  implement  such  a  technique  is 
shown  in  the  end  of  the  thesis. 


Figure  6.1.   How  to  Define  the  Cost  Function 

For  the  ideal  case  (Q=I  and  P  an  accurate  transformation) 
SQUERR  was  found  to  be  .0  966  and  for  the  non-ideal  one  worked 
in  run  No.  10  (Q^I  and  same  P  as  before)  it  was  computed  as 
.1117.   If  the  compensated  Fourier  coefficients  are  used 
instead  of  the  gradient  search  in  any  of  the  given  non-ideal 
cases,  the  cost  becomes  much  higher  (1.1  in  run  No.  10) .   In 
all  the  examples  worked  out  M  was  taken  as  twenty-one. 
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Another  important  feature  of  the  new  technique  is  that  it 

permits,  through  a  convenient  weighting  of  the  different  terms 

in  the  functional  cost,  the  desired  output  to  be  fitted  in 

certain  zones  even  more  exactly  than  in  the  ideal  case.   One 

example  of  this  situation  is  given  in  the  run  No.  9  for  which 

SQUERR  was  defined  as 

M 
SQUERR  =   J  k. (yl.  -  y2.) 
i=l  1   L     1 

where   M   is    as   before    and 

k.       =      1    for   i=l    to    7    and    15    to   2 1 

i 

k.   =10  for  i=8  to  10  and  12  to  14 

l 

k.   =   100  for  i=ll 

l 

1 .   Run  No.  8  Saturated  Amplifier  with   ~  ~ 

Up  to  now  the  signal  of  the  flux  given  by  the  manipu- 
lators was  not  restricted.   This  does  not  mean  that  the  flux 
may  be  negative  (removal  of  heat)  but  that  a  change  of  variable 
was  done  in  the  equation  describing  the  process,  such  that  the 
simulated  heat  flow  corresponds  physically  to  a  higher  level 
heat  flow.   In  some  circumstances  the  mentioned  transforma- 
tion may  be  undesirable  and  therefore  it  becomes  necessary  to 
take  into  account  in  the  simulation  the  fact  that  the  control 
cannot  go  negative.   The  output  was  far  from  the  desired  one 
and  then  it  was  thought  that  maybe  another  set  of  x ' s  would 
be  better  in  nonlinear  cases. 

To  verify  the  suppositions  just  mentioned  and  knowing 
that  the  Fourier  coefficients  give  a  minimum  square  error 
distribution,  a  cost  function  was  defined  (SQUERR)  which  was 
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WANT 


1       X 


a)       T   =    0.0025 


b)       T   =    0.005 


1      X 


WANT 


c)       T    =    0.1 


1       X 


Figure     6.2.        Run   No. 
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the  summation  of  the  squares  of  the  differences  between  the 
desired  and  the  obtained  distributions  at  every  one -twentieth 
fraction  of  the  system's  dimension.   The  new  x's  did  not  come 
out  very  different  but  the  deviation  suggested  that  in  some 
cases  it  is  possible  to  have  a  better  output  with  a  set  of 
t's  different  from  the  Fourier  coefficients.   This  is  true, 
as  it  will  be  seen  in  the  next  runs. 

The  optimal  set  of  reference  coefficients  is 

TAU(l)   =   114.62 

TAU(2)   =    84.38 

TAU(3)   =    38.70 
The  remaining  of  the  program  is  as  in  run  No.  5,  except  for 
the  nonlinearity  which  is  defined  as 

40  DO  41  1=1, MP 

IF  (ALPHA ( I) .GT. 40.0)  ALPHA ( I) =40 . 0 

41  IF  (ALPHA(I) .LT.0.0)   ALPHA(I)=  0.0 
2.   Run  No.  9:   A  Locally  Better  Output  Distribution 

The  nonlinearity  was  removed  and  it  was  decided  to 
reduce  the  deviation  between  the  theoretic  and  real  outputs 
in  the  central  zone.   In  order  to  achieve  this  result  a 
greater  weight  was  given  to  the  differences  between  positions 
eight  and  fourteen,  namely  ten  times  for  all  except  the  central 
one  (number  eleven)  which  received  a  weight  of  one  hundred.   As 
a  result  the  square  error  is  not  minimized  but  the  distribution 
approaches  the  desired  one. 
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WANT 


1      x 


a)       T  =    0.0025 


b)       T   =    0.005 


1       X 


c)       T    =    0.1 


1        X 


Figure     6.3.        Run   No.    9 
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The  resultant  set  of  t's  was: 
TAU(l)   =   116.9 
TAU(2)   =    88.2 
TAU(3)      -35.6 
The  consequences  of  these  results  may  be  quite  impor- 
tant.  Although  in  the  present  problem  three  ei gen functions 
approximate  already  closely  the  required  output,  in  other 
circumstances  this  may  not  happen.   Then,  the  possibility  of 
being  able  to  get  a  closer  output  in  some  regions  may  indeed 
be  relevant. 

3.   Run  No.  10:   A  different  Manipulator  Control 

The  computer  program  for  this  case  is  the  one  given  as 
reference.   As  mentioned  in  the  beginning  of  this  section  this 
is  a  situation  for  which  the  use  of  the  compensated  coeffi- 
cients is  not  recommended.   As  it  can  be  seen  in  the  program 
the  coefficients  obtained  using  the  gradient  search  came  out 
substantially  different  from  the  theoretical  ones. 
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WANT 


a)       T  =    0.0025 


1      X 


b)       T    =    0.005 


1       X 


HAVE 


c)       T    =    0.1 


Figure      6.4.        Run   No.     10 
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VII.       CONCLUDING    REMARKS 

A.       SUGGESTIONS    FOR   FUTURE    RESEARCH 

Many   more    difficult   problems    than   the    one   worked   out   may 
appear,    such   as    cases    involving: 

i)       Non-homogeneous    boundary    conditions;    a  transformation 
of   variable    is    very    useful   in    cases    like    these. 
ii)       Discontinuities, 
iii)       Time    varying    coefficients, 
iv)       Hyperbolic   differential   equations.       In   this    case,    if 
using   the   Bubnov-Galerkin   transformation   the    set    of 
state    variables   must    be    augmented  with   those    correspon- 
dent to    the    second   derivative.       Such   procedure    implies 
the   need    for  knowing   the    initial    conditions    of   the    first 
derivatives,    which   is   not    a   normal    situation,    and   there- 
fore  the    use    of    an   observer   has    to   be    implemented. 
For   the    beginner   that   gets    involved   in   difficult    situa- 
tions   it  may   be   necessary   to    recur   to    the    applied   mathemati- 
cian.     However,    before    that,    there    are    three    references    that 
must   be    considered   very    carefully   because    of   the   tremendous 
potential   of   knowledge    they    contain:       Courant    and  Hilbert 
[11],    Crandall    [12]    and  Murray-Lasso    [45].      Also   Foster    [18], 
working   in    an    optimal   control    fashion,    develops    a   systematic 
procedure    for   the   modeling   of    linear   regulator  problems, 
using   the   Bubnov-Galerkin    transformation. 
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As  areas  for  future  research,  problems  involving  non- 
linearities  in  the  control,  time  variant  coefficients  and 
parameter  identification  are  recommended,  as  well  as  a  gener- 
alization of  the  Foster  optimal  control  method  for  a  class  of 
problems  other  than  the  linear  regulator. 

B.   CONCLUSIONS 

The  most  significant  of  the  original  contributions  con- 
tained in  the  present  work  is,  certainly,  the  first  closed- 
loop  simulation  of  the  modal  control  of  a  distributed 
parameter  system  by  non-optimal  control  techniques.   The  con- 
clusions may  be  summarized  as  follows: 

i)   The  classical  control  methods  can  be  applied  to  each 

independent  loop.   A  constant  in  the  feedback  is  general- 
ly good  enough  to  obtain  faster  response.   If  it  is 
desired  to  work  close  to  the  stability  limit,  design 
procedures  such  as  Bode  diagrams  or  root-locus  plots 
also  seem  to  be  very  suitable  to  this  purpose. 
ii)   The  available  theory  states  that  whenever  the  control  is 
linear  (at  least  in  steady-state) ,  the  Fourier  Coeffi- 
cients of  the  desired  output  distribution,  multiplied 

by  the  factor  -^ — ^— - — *■  ,  give  the  compensated  Fourier 

D  3 

coefficients  which  used  as  forcing  functions  (t)  insure 

that  the  mean  square  error  of  the  output  is  minimized. 
As  shown  in  the  computer  results  the  above  state- 
ment was  proven  to  be  not  always  true,  specifically  when 
the  transformations  Q  or  P  are  not  accurate,  which  is  the 
most  general  situation.   Notice  however  that  Q  is  an 
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exact  transformation  (identity  matrix)  for  the  case  of 
eigenfunction  control.   Also  P  and  Q  may  be  quite  inac- 
curate when  they  are  low  order  matrices  and  it  is  neces- 
sary to  use  a  large  number  of  eigenfunctions  to 
approximate  with  precision  the  output  or  the  control. 
When  this  happens  the  compensated  coefficients  may  be 
quite  far  from  optimal;  such  situation  can  be  partially 
avoided  using  the  technique  of  gradient  search  for  the 
reference  coefficients  x,  as  originally  developed  in  this 
thesis, 
iii)   Through  an  adequate  definition  of  a  cost  function  it  is 
possible  to  vary  the  reference  coefficients  in  such  a 
way  that  the  minimum  square  error  is  always  achieved  or 
even  such  that  some  regions  fit  the  required  distribu- 
tion more  closely  than  that  minimum  square  error 
approximation . 
iv)   The  best  set  of  manipulators  is  the  one  that  gives  a 
heat  distribution  coincident  with  the  eigenfunctions. 
In  this  case  G  becomes  the  identity  matrix  and  the  sys- 
tem can  be  easily  modeled  using  the  Bubnov-Galerkin 
method.   In  many  circumstances  it  will  be  difficult  to 
have  that  type  of  manipulators;  then,  the  set  that  ap- 
proximates more  closely  that  distribution  is  the  desired 
one. 
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APPENDIX    A 

TRANSFORMATION    OF    THE    PARTIAL    DIFFERENTIAL    EQUATION    INTO    AN 
ORDINARY    DIFFERENTIAL    EQUATION    BY    THE    BUBNOV-GALERKIN    METHOD 

The   Bubnov-Galerkin   method   much    used   by   Foster    [Ref .    17] , 
and   described   in   great    detail   by   Mikhlin    and    Smolitskiy    [Ref. 
43] ,    is    an    approximation   technique    for   the    solution    of   the 
equation  N   y-q=0,    where   N    is    a   differential    operator    (•—  -   L 

o  t  X 

in  this  case)  that  must  obey  certain  conditions.   The  condi- 
tions are  boundedness  and  existence  of  a  completely  continuous 
inverse,  and  they  will  be  satisfied  if  y(x,t)  is  unique  and 
the  coefficients  of  N,  as  well  as  their  first  derivatives  are 
continuous . 

a 

Although  apparently  limited,  the  technique  applies  to  a 
large  variety  of  physical  situations  and  it  seems  to  start 
playing  a  very  important  role  in  the  control  of  D.P.S.'s.   It 
is  a  generalized  and  systematic  way  of  getting  the  uncoupling 
effect  obtained  in  Chapter  IV. 

Start  by  writing  an  approximate  solution  of  the  partial 

differential  equation  in  terms  of  the  eigenfunctions : 

N  18 

y(x,t)    =      I    co    (t)    u    (x)     =      I    oj    (t)    un(x)  (A.l) 

n=l  n=l 


1  o 

In   the    simulation   run  No.    7    vector   oo   is    represented 
by   vector   Z . 
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Applying  the  condition  that  Ny  -  q  is  orthogonal  to  the 
N  eigenfunctions  considered  here  it  follows  that 

do>(t) 


dt 


=  Au)(t)  +  £(t)  (A. 2) 


where 


oj(t)  =  [Ul(t)  o)2(t)  —  (^(t)]7 
C(t)  =  [^(t)  c2(t)  —  cN(t)]T 

C±(t)  =  <q(x,t)  ,  ui(x)> 


[A]  .  .  =  <u.  (x)  ,  L   u.  (x)  > 

J-  j       i        x   j 


1.   Solving  for  £. (t) 


but, 


(A.  3) 


1 
5±(t)  =  /  q(x,t)  u±(x)  dx  (A. 4) 


M 
q(x,t)  =  I    am(t)  Hm  (x,£m)  (A. 5) 

m=l 

where  a   and  H   were  defined  in  Eqs .  4.56  and  4.51.   Substi- 

m      m 

tuting  A. 5  in  A. 4  it  follows 

M 

?i(t)  -  I.    bmi  «m(t>  (A-6) 

m=l 

with 

1 

b.  =  /  u,  (x)  H   (x,  e  )  dx  (A. 7) 

mi    A 

which  is  the  same  as  the  previous  definition  of  the  elements 
of  B. 

Therefore  Eq.  A. 2  may  now  be  rewritten  as 

du(t)  T 

— ^r =  Aco(t)  +  B  a  (t)  (A.  8) 

at      ~~       ~ 
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In  the  case  of  N  filters  and  manipulators  ct(t)  = 

(B  )   p.  and  B  a(t)  =  BT  (BT)  ~\    =  y  which  means  that  it  is 

T 
not  necessary  to  consider  the  matrix  B   in  the  model.   As  a 

consequence,  the  simulation  by  this  method  only  reproduces 
the  real  situation  when  the  heat  distribution  of  the  manip- 
ulators coincides  with  the  eigenfunctions,  giving  B  =  I. 
2.   Solving  for  [A] .  . 


[A]    =  /  u.  (x)  -^4  dx  (A. 9) 

1D    0   1       8xZ 

and,  from  Chapter  V,  u.  (x,s)  =  /2    sin  A-s  x.   Because  M  and 

+ 
M   are  self -ad joints ,  the  normalized  eigenfunction  set  is 

orthogonal  and  applying  the  inherent  properties  the  matrix  A 
becomes  diagonal. 

A  parenthesis  is  opened  to  say  that  the  computation  of 
the  eigenfunctions  done  in  Chapter  V.B  was  not  indispensable. 
As  a  matter  of  fact  it  would  have  been  enough  to  choose  an 
arbitrary  orthogonal  basis  forming  a  complete  set  and  satisfy- 
ing the  boundary  conditions,  and  do  all  the  derivations  with 
this  basis.   One  set  that  is  chosen  very  often  and  that,  by 
coincidence,  is  the  derived  set  of  eigenfunctions  is 

3±(x)  =  /2  sin  ^ix  .  .  (A.  10) 

Back  to  the  elements  of  A  their  values  are: 
A±1      =   -  |_   Sin2  |x  dx  _  _  n_  =  _2.46740 

9J12    r1   .    2   3Hx  ,      9n2    00  „,,,, 
A?9   =   ~  0 /  Sln  ~T~   dx  =  "  4 =  -22.20661 


*22 


0  (A. 11) 


25II2      r1     .      2     5ITX      ,  25H2  cl      CQrn-5 

A.....      =      -   —^ J    sin      —jT—  dx   =    -   — j =    -    61.68503 


33 

A.  .       =    0,    i    i-    j 
13  1         7     j 


0 
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and  now  all  the  elements  in  Eq.  A. 8  are  known.   The  initial 
conditions  are  £■  (0)  ,  equal  to  the  i    coefficient  of  the  ex- 
pansion of  the  initial  state  y(xr0)  and  for  the  specific 
problem  treated  here  they  turn  out  to  be  zero. 
The  output  is  given  by 

y(x,t)  =  C(x)u)(t)  (A. 12) 

where 

C(x)  =  [u±(x)    u2(x)  ...  i^Cx)]        (A. 13) 
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APPENDIX  B 
BASIC  DEFINITIONS  IN  FUNCTIONAL  ANALYSIS 


1.   Operator 

An  operator  is  a  mapping  of  functions  into  functions. 
Figure  B.l  is  a  block  diagram  representation  of  the  operator 
L,  which  maps  the  function  m(x,t)  into  the  function  y(x,t). 
The  set  {m}  for  which  L  is  defined  is  the  domain  of  L.   The 
corresponding  set  {y}  is  the  range  of  L.   An  operator  can 
take  several  configurations,  the  most  frequent  of  which  are 
the  differential  operator  and  its  inverse,  the  integral 
operator. 

When  computing  the  inverse  of  an  operator  (m(x,t)  = 
L  y(x,t))  it  is  necessary  to  know  the  points  at  which  it  is 
singular  and  these  points  are  called  the  eigenvalues  X.,    of 
the  direct  operator  (L) .   The  eigenf unctions  are  then  defined 
as  the  functions  u  which  satisfy  the  equality 


L~u   =   Mu   =   X  u 


(B.l) 


m(x/t) 


Figure  B.l 


Representation  of  the  Uncontrolled  System 
in  Operator  Form 
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The  partial  differential  equations  considered  here  have 
an  infinite  complete  set  of  eigenf unctions  and  by  complete, 
it  is  meant  that  any  function  with  a  finite  number  of  dis- 
continuities can  be  represented  by  a  linear  combination  of 
these  eigenfunctions .   This  is  the  concept  of  generalized 
Fourier  series. 
2 .   Inner  Product 

The  inner  product  of  two  continuous  infinite  vectors  p 

and  q  is  defined  as 

b 
<p,q>=/pqdx  (B.2) 

a 

for  finite  values  of  the  integral.   If  the  vectors  turn  out 

to  be  discrete  infinite  this  becomes 

<P^>  =  N  +  oo    I     P^i  (B'3) 

r=l 
provided  that  also  this  summation  is  finite.   The  limit  is 
dropped  for  the  case  of  finite  vectors. 

When  two  sets  of  fectors  are  such  that 

<p,q>  =  0,  p  J   q  (B.4) 

=  1,  p  =  q 
these  two  sets  are  said  to  be  orthonormal. 

3 .   Adjoint  Operator 

+ 

The  adjoint  of  an  operator  M  is  the  operator  M   such  that 

<M  u,  v>  =  <u,  M  v>  (B.5) 

An  important  property  of  the  adjoint  operator  is  that  its 

eigenfunctions  (v.)  form  a  set  orthogonal  to  the  eigenfunc- 
tion's  set  (u. )  of  M. 
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If  M  =  M   the  operator  is  said  to  be  self-adjoint  and  in 
this  case  both  M  and  M  have  the  same  orthogonal  set  of 
eigenf unctions .   Its  properties  are  analogous  to  the  proper- 
ties of  the  symmetric  matrix  and  it  can  be  proved  that  the 
self-adjoint  case  gives  real  eigenvalues. 
4 .   Spectral  Representation  of  an  Operator 

Using  the  orthogonality  of  the  sets  u.  and  v. ,  previously 
normalized,  it  follows  easily  that  any  operator  with  a  dis- 
crete infinite  spectrum  may  be  represented  as 

b   °° 
M-]   =   /  {  JX.  u.  (x)vi(c)}[-]dc  (B.6) 

a   i=l 

which  takes  the  name  of  spectral  representation  of  an 

operator. 

If  q(x)  =  Lp(x),  from  the  last  equation  and  also  from  the 

OO  00 

Fourier  series  for   p(x)  =   £  a.  u. (x)  and  q(x)  =   £  b.  u. (x) 

i=l     L  i=l 

it  turns  out  that 


q(x)  =  I    a   L[u •  (x)]  =   J  X  .  a   u  (x) 
i=l  1     L       i=l 


(B.7) 


Therefore 


(B.8) 


I   b   u  (x)  =  I   X±a  u  (x) 
i=l  x      x  i=l 

from  which,  by  the  orthonormality  properties, 

b.   =   X.a.  (B.9) 

1  D  D 

This  means  that  the  coefficients  of  the  eigenfunction  ex- 
pansion of  the  output  can  simply  be  obtained  from  the  corres- 
ponding coefficients  of  the  input  by  a  simple  product  and 
therefore  the  operator  was  diagonalized. 
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A  useful   property   of    the    diagonalized   operator   is    the 

function   of   an   operator   which   simply   states 

b        °° 
f(L[.])       =      /     (      I    f(Xi)u.(x)v.(C))  [-]dC  (B.10) 

a        i=l 

For  the  important  case  of  the  inverse  operator  the  spectral 

representation  is 

b  °° 
!•   !•]   =  I'il      r~ u,  (x)v.  (C)>[-]dC  (B.ll) 

a   i=l   i  1    x 

In  Ref.  45  and  extensive  theory  of  operators  is  given , 
which  furnishes  a  deep  knowledge  about  procedures  for  dealing 
with  complex  situations.   This  thesis  will  be  restricted  to 
separable  operators  for  which  in  the  case  of  partial  differen- 
tial operators  of  the  type 

L   . Im(x,t)]   =   H. [m]  +  H  [m]  (B.12] 

x  i  l.  u-         a. 

The    separability    is    realized   if 

Ht[Hx[-]]       =      Hx    [Ht[-]]  (B.13) 

Further    restricting  H      to  have    constant   coefficients    and 
Laplace   transforming   the   time,    the    spectral    representation 
of    L  becomes 

1 

L         !•]       =         I    A.  (s)u.  (x)J       v.      (C)     dC  (B.14) 

x,s  i=1    x  x         0         x 

A  more  systematic  way  of  obtaining  the  diagonalization  of 
an  operator  is  using  the  Bubnov-Galerkin  transformation,  as 
described  in  Appendix  A. 


146 


5 .   Other  Definitions 

Linear  operator:   an  operator  L  is  linear  if  the  mapping 
that  it  implies  is  such  that  for  arbitrary  scalars  a  and  b 
it  follows  that  L(ax,  +  bx«)  =  aLx-j+bLx^ . 

Continuous  operator:  a  continuous  operator  is  character- 
ized by  the  fact  that  if  a  sequence  of  vectors  {x  }  converges 
to  x,  then  the  sequence  of  vectors  {Lx  }  converges  to  Lx. 

Completely  continuous  operator:   a  linear  operator  is 
completely  continuous  if  for  every  bounded  sequence  {x  }  in 
a  linear  normed  space  X  the  set  {Lx>   has  a  subsequence  which 
converges  in  the  mean  to  an  element  of  X. 

Convergence  in  the  Mean,  also  called  strong  convergence: 
the  sequence  {x  }  is  said  to  converge  in  the  mean  to  x  when 


2 

L   space,  which  implies  that  x  is  a  function  defined  in  the 


interval  a-b)  or 


/      2 
x   =  \      1    |x.  |<oo  (for  {x}  a  set  of  complex 

V  i=l   1 


2 

numbers  and  it  is  said  that  the  space  is  an  1   space) . 

Bounded  operator:   a  linear  operator  L  is  bounded  if  its 
domain  is  the  entire  space  and  if  there  is  a  scalar  M  such 


that 


Lxi  <  M  li  x  I  •   The  smallest  of  these  bounds  is  called 


norm  of  the  operator  and  denoted  by  i;  M 

Subsequence :   given  a  sequence  S:   {Lx,  ,  Lx~  ,  Lx^  ,...} 

a  subsequence  of  S  is  a  fraction  of  it  given  by  Lx   ,  Lx   , 

nl     n2 

•3 


Lx  ,    ...Lx    ...   where  n, <n~ <n_< . . .n . < . . .  . 
n-3        n^  1   2   3      l 
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c 
c 

c 
c 

C  BROAD  COMPUTATION  OF  THE  OPTIMAL  POSITIONS  FOR 

C      -   THE  SENSORS. 

C  THIS  PROGRAM  DOESN'T  TAKE  INTO  ACCOUNT  THE 

C  POSSIBILITY  OF  ONE  SENSOR  AT  X=0.0  BECAUSE 

C  THE  TEMPERATURE  IS  KNOWN  AT  THIS  POINT  FROM 

C  THE  BOUNDARY  CONDITIONS. 

C 

C 

C 

C  THE  MAIN  PROGRAM  MUST  BE  ADJUSTED  ACCORDING  TO  THE 

C  NUMBER  OF  SENSORS. HOWEVER , THE  SUBPROGRAM  USED 

C  (EVMAX)  IS  VERY  GENERAL  AND  REQUIRES  ONLY  THE 

C  ADJUSTMENT  OF  THE  DIMENSIONS  AND  DATA  AND  OF 

C  THE  LOOP  "DO  10  ». 

C  EVMAX    DOES    NOT    HOLD    IF    MATRIX    A    IS    COMPLEX 

C 

C 

c 

DIMENSION    X(3) ,Y13) 

EXTERNAL  EVMAX 
C 

C      N=   NUMBER  OF  MEASUREMENTS 
C      NG=   NUMBER  OF  DIVISIONS 
C 
C 

DATA  N/3/,NG/10/ 

Z=.1E25 

1  NGl=NG-2 
NG2=NG-1 
NG3=NG 
Nl  =  l 
N2=N1 

DO  2  I1=N1,NG1 

N2=N2+1 

N3  =  N2 

DO  2  I2=N2,NG2 

N3=N3+1 

DO  2  I3=N3,NG3 

X  (  1 )  =  (  .  1  )  *  1 1 

X(2)=( .1)*I2 

X(3)=(  .1)*I3 

IF(EVMAX(X) .LT.ZJGO  TO  3 

GO  TO  2 

3  DO  4  I  =  1,N 

4  Y(I)=X( I  ) 
Z=EVMAX(X) 

2  CONTINUE 

7  WRITE  (6,5) 
C 

C      NEXT  THREE  STATEMENTS  MUST  BE  ADJUSTED  FOR  N.GT.3 
C 

WRITE  (6,6)  Z, ( Y( I ) ,1=1, N) 

5  FORMAT  ( '1« ,////////, T35, 'EVMAX* ,T49, 'XI1 ,T64,« X2» , 
1T79, 'X3' ,////) 

6  FORMAT  ('  « ,T33,E10.3, T42, 3 ( F 10 . 3 , 5X ) ) 
STOP 

END 

FUNCTION  EVMAX(X) 
C 
C 

C  THE  DIMENSION  MUST  BE  ADJUSTED  IN  THE  NEXT 

C  STATEMENT  ACCORDING  TO  THE  NUMBER  OF  EIGENVALUES 

C  THAT  APPROXIMATE  THE  SYSTEM. 

C  ALSO  THE  DATA  MUST  BE  CHANGED  IF  A  DIFFERENT 

C  SYSTEM  IS  USED. 

C  EVMAX  DOESN'T  HOLD  IF  MATRIX  "A"  IS  COMPLEX. 
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c 
c 


c 
c 
c 
c 
c 
c 

c 
c 
c 
c 


10 
15 


30 


DIMENSION  A(3),X(3),U(3,3),UI(3,3),L(3),M(3),UIT(3,3), 
1EVABS(3) ,UIUIT(3,3),EIVU(3) 
DATA  N/3/, A/ 1.57 08 0,4. 7123  9, 7.85398/ 
DO  10  J=1,N 
DO  10  1=1, N 
B=A( J)*X(  I  ) 
U(  I,J)=1.41421*SIN(B) 
UK  I,J)=U(  I, J) 
DO  15  1=1, N 
WRITE(6,300) 
CALL  MINVtUI 
IF(D.EQ.0.0)GC 
DO  30  J=1,N 
DO  30  1=1, N 
UIT( J, I )=UI ( I, 


20 


25 


40 


46 
45 

50 
200 
300 
400 

60 


I 
»N 


,X(I  ) 

,D,L,M) 


TO 


3  ' 


J) 


CALL  MPRD(UIT,UI,UIUIT,N,N,0,0,N) 

DO  20  1=1, N 

DO  20  J=1,N 

WRITE  (6,200)  I , J,UIUIT(I, J) 

S/R  JACVAT  COMPUTES  THE  EIGENVALUES  OF  A  PEAL 
SYMMETRIC  MATRIX, WHICH  IS  OUR  CASE, BECAUSE  UIUIT 
IS  THE  PRODUCT  OF   MATRIX  "UIT"  AND  ITS  TRANSPOSE 
"UI" 

CALL  JACVAT  ( J IU IT , N, 0, E I VU, DUMMY, N ) 

EIVU  IS  THE  VECTOR  OF  COMPUTED  EIGENVALUES  AND 
EVMAX  IS  THE  MAXIMUM  EIGENVALUE 

DO  25  1=1, N 

WRITE  (6,400)  I,EIVU(I) 

DO  40  1=1, N 

EVABSd  )=ABS(EIVU(  I)  ) 

CONTINUE 

EVMAX=-1.0 

DO  45  1=1, N 

IF(EVABS(I ) .GT. EVMAX)  GO  TO  46 

EVMAX= EVMAX 

GO   TO  45 

EVMAX=EVABS(I) 

CONTINUE 

GO  TO  60 

EVMAX=. 1E2  0 


FORMAT 
FORMAT 
FORMAT 
RETURN 
END 


•0' 
•0' 
'0' 


"JIUIT( 
•  X  (  »  ,  I  2 
•EIVU( • 


t  I 
•  ) 
12, 


2, 
—  i 


t',  12, 
F10.5) 


•  )=' ,D15.5) 


•,E14.5) 
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c 

c 
c 
c 

C         RIGOROUS  COMPUTATION  OF  THE  SENSORS   

C         POSITIONS  BY  GRADIENT  SEARCH  

C 

c 
c 
c 

C         THE  MAIN  PROGRAM  AND  THE  SUBPROGRAM  EVMAX 

C         ARE  VERY  GENERAL  AND  THE  ONLY  THING  THAT  NEEDS 

C         TO  BE  CHANGED  ACCORDING  TO  THE  SYSTEM  ARE  THE 

C         DIMENSIONS, THE  DATA, LOOP  "DO  10 "  AND, IF 

C         NECESSARY, THE  PARAMETERS  OF  S/R  DIRECT. 

C         EVMAX  DOESN'T  HOLD  IF  MATRIX  "A"  IS  COMPLEX. 

C 

C 

DIMENSION  X(3) 
EXTERNAL  EVMAX 
DATA  N/3/,X/.3, .6, .9/ 
DELCAP=.01 

CALL  DIRECT  ( X ,N, EVMIN, DELCAP ,. 25, 1 . E-4, EVMAX, KONVRG , 
10,1) 
WRITE  (6,300) 
WRITE  (6,400)  KONVRG 
WRITE  (6,500) 
DO  10  1=1, N 
10  WRITE  (6,600)  X( I) 
300  FORMAT ( • 1' ,T60, 'KONVRG' ) 
400  FORMAT  ( • 0 • , T60, 14 , //// ) 
500  FORMAT  ( ■ 0 ' ,T60 , • X ' ) 
600  FORMAT  ( • 0 ' ,T50 , Fl 5 .5 ) 
STOP 
END 

FUNCTION  EVMAX(X) 
C 
C 

C         THE  DIMENSION  MUST  BE  ADJUSTED  IN  THE  NEXT  THREE 
C         CARDS  ACCORDING  TO  THE  NUMBER  OF  EIGENVALUES 
C         THAT  APPROXIMATE  THE  SYSTEM. 

C         ALSO  THE  DATA  MUST  BE  CHANGED  IF  A  DIFFERENT 
C         SYSTEM  IS  USED. 
C 
C 

DIMENSION  A(3)  ,X(3)  ,U(3,3) ,UI( 3,3) ,L(3) ,M(3) ,UIT(3,3) , 
1EVABS( 3),UIUIT(3,3 ) ,DI VU13) 

DATA  N/3/, A/ 1.5708 0,4. 71239, 7. 85398/ 

DO  10  J=1,N 

DO  10  1=1, N 

6=A(J)*X(I ) 

U(  I  ,J)=1.41421*SIN(B) 
10  UK  I,  J)  =  U(  I  ,J) 

DO  15  1=1, N 
15  WRITE(o,300)  I ,X(I  ) 

CALL  MINV( UI ,N,P,L ,M) 

IF(D.EQ.0.0)G0  TO  50 

DO  30  J=1,N 

DO  30  1=1, N 
30  UIT(J, I )=UI( I,  J) 

CALL  MPRD(UIT, UI , UI UIT ,N,N ,0 , 0 , N ) 

DO  20  1=1, N 

DO  20  J=1,N 
20  WRITE  (6,200)  I , J, U IUI T( I , J ) 
C 

C         S/R  JACVAT  COMPUTES  THE  EIGENVALUES  OF  A  REAL 
C         SYMMETRIC  MATRIX, WHICH  IS  OUR  CASE, BECAUSE  UIUIT 
C  IS  THE  PRODUCT  OF   MATRIX  "UIT"  AND  ITS  TRANSPOSE 

C         "UI" 
C 
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CALL  JACVAT  ( U I UIT ,N, 0 , EI VU, DUMMY, N) 
C 

C      EIVU  IS  THE  VECTOR  OF  COMPUTED  EIGENVALUES  AND 
C      EVMAX  IS  THE  MAXIMUM  EIGENVALUE 


C 


-DO  2  5  I=1,N 
25  WRITE  (6,400)  I,EIVU(I) 
DO  40  1=1, N 
EVABSU  )=ABS(EIVU(  I)) 
40  CONTINUE 
EVMAX=-1.0 
DO  45  1=1, N 

IF(EVABS( I ).GT. EVMAX)  GO  TO  46 
EVMAX= EVMAX 
GO   TO  45 
46  EVMAX=  EVABSU) 
45  CONTINUE 
GO  TO  60 
50  EVMAX=.1E20 
200  FORMAT  ( • 0 • , ' U IU IT ( • , I  2 , ' , ' ,  12 , • )= » , D15 . 5) 
300  FORMAT  (  ' 0' , ' X ( '  , 12  ,  «  ) =  ■  , F 10 . 5 ) 
400  FORMAT  ( • 0 ■ , ' E I VU( ' , 12 , ■ )= ■ , E14. 5 ) 
60  RETURN 
END 
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c 
c 
c 

C  MODAL  CONTROL  OF  A  DISTRIBUTED  PARAMETER  SYSTEM 

C  DESCRIBED   BY  A  PARABOLIC  DIFFERENTIAL  EQUATION. 

C 

C  CARDS  UP  to  "no  4 — "  AND  THE  DATA  CARDS  MUST  RE 

C         CHANGED  ACCORDING  TO  WHAT  IS  REQUIRED  FROM  THE 

C  SYSTEM  AND  ALSO  ACCORDING  TO  THE  PHYSICAL  LIMITA 

C  TIONS.THE  3NLY  ADDITIONAL  CHANGES  MAY  3E  IN  LOOP' 

C         "DO  3« — ", WHICH  DEFINES  TMF  HEAT  DISTRIBUTION  HF 

C         THE  MANIPULATORS, IN  STATEMENT  #43,  IN  THE  LOOP 

C         "DO  41 — »  AND  IN  THE  PARAMETERS  AND  VECTOR  Y2  OF 

C         S/R  "DRAW". 

C 

C 

REAL*8  R,A1,ITITL1(12),ITITL2(12),ITITL3(12) 

REAL*4  MIU 

REAL  LABEL1/4HWANT/ 

REAL  LABEL2/4HHAVE/ 

DIMENSION  ALPHA(3) ,U(21,40   ),R(20),A(59  ),Q(3,3), 
1P(3,3),X(3),EV1(21),EV2(21),EV3(21),H1(21),H2(21), 
2H3(?1)  ,  TAU(3)  ,TEMPCM  ,X1(21)  ,  YH21)  ,°HI  (3)  ,  Al(5°)  , 
3RK80)  ,  R2(2C)  ,MIU(3)  ,0'1EGA(3)  ,  RO(  3  )  ,  Y2  (  2  1  ) 
C 

C      N #    OF  FQUATIONS 

c      m a    OF  RIGHT  HAND  VECTORS  IN  THE  SYSTEM  OF  EQUA_ 

C      TICNS  TO  BE  SOLVED  DY  CRANK-NI COLSON  METHOD. 

C      KF V-    OF  «=ILTFDS 

C      KS H    OF  SENSORS 

C      MP #  OF  MANIPULATORS 

C 

OATA  N,M,KF,KS,MP,NT/20,4,3,3,3,40   / 

XU)  =  . 28570 

X<2)=.  5_7]44 

X(3l=.85715 

TAU(1)=102  .86 

TAU( 2) -62. 15 

TAU(?) =-60.61 

PHI (] )=118.4353 

PHI (2) =°B.6961 

PHI (3)=5Q.2177 
C 

C      READING  THE  MATRICES  P  AND  Q  OF  THE  FEEDBACK  LOOP 
C 

DO  4  I=1,KF 

4  READ  (5,2000)  ( P ( I , J ) , J=l , KS) 
DO  5  1=1, MP 

5  READ  (5,2000)  (  Q  (  I  ,  J  )  ,  J  =  l  ,  !<F  ) 
N1=N+1 

N2=N+2 

H=1.0/N 

HK=H**2 

Jl  =  l 
C 

C      BOUNDARY  CONDITIONS  AT  X=0 
C 

DO  2  J=1,NT 

2  U  (  1 ,  J  )  =  0 .  0 
C 

C      INITIAL  CONDITIONS 
C 

DO  3  1=1, Nl 

3  UCItl )=0.0 
C 

C      SOLUTION  BY  CR ANK-NTCHQLSON  METHOD 
C 

c 

C      DEFINING  THE  TRIDIAGONAL  MATRIX  OF  EQUATIONS  ACCOR- 

C      DING  TO  THE  REQUIREMENTS  OF  S/R  DGELB. 

C 

NM=N*(M-1 ) 
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c 
c 

c 


c 
c 
c 
c 
c 
c 


c 
c 
c 
c 


c 
c 
c 
c 


c 
c 

c 


30 


NM0=NM-2 

N0=NM-3 

DO  30  1=1, NMO, 3 

A( I )=4.0 

IF(I.EQ.NMO)  30  TO  30 

11=1+1 

12=1+2 

A( Il)=-1.0 

A( I2)=-1.0 

CONTINUE 

A(MC)=-2 


0 

DEFINING  THE  FIRST  (M-1)*N  TERMS  0^  Rl 

DO  3?  IA=1,NM 
35  R1(IA)=0.0 

HEAT  DISTRIBUTION  Oc  THE  MAN  I PUL ATOP S . THF  NEXT  CARDS, 
UP  TO  STATEMENT  *38  MUST  BE  CHANGED  IF  THE  DISTRIBU- 
TION OF  HEAT  CORRESPONDENT  TO  EACH  MANIPULATOR  ALSO 
CHANGES. 


DO    38    ID=1,N1 


1) 


38 


EVK  ID)  =  6.283186-J-'H*(  ID 
IF(  IP.GF.1DFVK  ID)=P. 
EV?( ID)=3.141593*H*(ID-1) 
EV3(  ID)=3.1415<>3*H*<  ID-1  ) 
IF     (  I D  .  L E  .  1 1 )     E V 3  (  I  D )  =  0 . 
Hl( ID)=STN(EV1 ( ID) ) 
H?(ID)=SIN(EV2 (ID) ) 
H3( ID)=-C0S(FV3( IP) ) 

IF    I.C.'S.NE.O     INSTEAD    OF    TAU    WE    MUST    COMPUTE    AMD    USE 
MIU. 


CALL    MDRP(C\TAU,  AL PHA , MP , KF , 0 , 0,  1  ) 
WPITE     (6,1000)     ALPHA 
40    CONTINUE 

42  DO    43    1=1  ,N 
10=1+1 
IE=NM+I 

43  Rl(  IF )  =  ( ALPHA  (  1)*H1 (  ID ) +AL PHA ( 2 ) *H2 (  I D ) +AL PHA ( 3 ) * 
1H3UD1  )*HK*2.0 

WPITE(6,1000)     ALPHA 

S/R  DGELB  SOLVES  IN  DOUBLE  PRECISION  A  SYSTEM  OF 
LINEAR  EQUATIONS. 

DO  44  1=1, MMO 

44  Aid  )=A(  I) 
N3=2*N+1 
N4=3*N+1 

DO  4  5  1=1, N 

R2( I )=R1 ( I )+Rl ( Nl ) +R 1 ( N3 ) +R1 ( N4 ) 

M1=N1+1 
N3=N3+I 

45  N4=N4+1 

DO  47  I=1,N 
R(I )=P?(  I) 

47  CONTINUE 
MUD=1 
MLD=1 

CALL  DGELB  { R , A  1 ,N , 1 , MUD ,MLD, . 1E-06, I FR1 ) 
WPITE(6,8000)  IFR1 

REDEFINING  THE  I  .C.  'S 

DO  48  ID=2,M 

48  R1(ID)=0.0 
NE=2*N-1 
N0=N-1 

NE1=ME+1 
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c 
c 
c 
c 


IG=1 

DO    49    ID=N2,NE 
R1(ID)=R(IGJ 
49.    IG=IG  +  1 

-R1(NE1)=2*R(N0) 

NF=3*N 

NE2=NF+2 

IG  =  2 

DO    150     ID  =  NE2,\|F 

R1(ID)=R(IG) 
50  IG=TG+1 

R1(NF)=0.0 

N1=N+1 

DO  60  IP=2,N1 

TF=in-i 

U(TD,J1)=R(IE) 
60  WRITE  (6,3000)  I D, J  1 , U ( ID, J  1  ) 

COMPUTATION  OF  VECTORS  ARG  AND  VAL  FOR  THE  INTERPO_ 
LATING  S/R  "ALI". 

CALL  ANA  (TEMP,KS, J  1 ,H , N ,U , X , Nl , NT ) 

CALL  FEED  (  TEM P , AL PHA , KF , KS , MP , P ,Q , PHI , TA U ,MIU, OMEGA, 
1R0) 
WRITE  (6,6000)  J  1 , ( AL PHA ( I )  , I  =  1 , MP) 
J1=J1+1 

IF(Jl.GT.NT)  GO  TO  90 
GO  TO  4  0 

90  DO  °1  TH=1,N1 
Xl( IH)=( IH-1 )*H 
D=6.283185*X1( IH) 
Y2(  IH)  =  1.0-CPS (D) 

91  CONTINUE 
10=1 

DP  92  IH=1,N1 

92  Yl( IH)=U(TH, 13) 
PEAD(5,4000)  ITITL1 
PEAP(5,4000)  ITITL2 
PEAP(5,4000)  ITITL3 

CALL  DK AW( 2 1, XI, Y2, 1,0, LABEL  1, I TITL 1,0,0, 1,1,2,2,7,7, 
11, LAST) 

CALL  DRAW(21,X1,Y1,3,0,LABEL2,ITITL1,0,0,1,1,2,2,7,7, 
11, LAST) 

IQ  =  2 

DO  93  IH=1,N1 

93  Yl( IH)=U(IH, 10) 

CALL  DPAW(21 , X 1 , Y2 , 1 , 0, LABEL  1 , I T I TL2 ,0 , 0 , 1 , 1 , 2 , 2 , 7, 7, 
11, LAST) 

CALL  DPAWm  ,X1,Y1,3,0,LABEL2,ITITL2,0,0,  1,  1,2,2,7,7, 
11, LAST) 

DO  94  IH=1,N1 

10  =  40 

94  Yl( IH)=U( IH, 10) 

CALL  DP.  AW(  2  1,X1,Y2,  1,0,  LAB  ELI,  I TITL3,  0,0,  1,1,2,2,7,7, 
11, LAST) 

CALL  DRAW( 21,X1, Yl , 3, 0, LABEL  2 , I TITL3 ,0, 0 , 1 , 1 , 2 ,2 ,7 , 7 , 
'  1,LAST) 

FORMAT  ( «0' , • ALPHA=' ,3E15.5) 

FORMAT  (8E10.5) 


1000 
2000 
3000 
4000 
6000 
8000 


) 


NE15.5) 


C 

c 


FORMAT  ( • 

FOPMAT( 10A8) 

FORMAT  (»0','ALPHA(,,I6,,)=,,7E14.5) 

FORMAT  (  '0'  ,  MEP1='  ,15  ) 

END 
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SUBROUTINE    ANA    ( TEMP* KS, Jl , H,N, Ut X,N1 »NT) 
C 

c 

C      THIS  S/R  COMPUTES  THE  INTERPOLATED  VALUES  OF  THE 
C      TEMPERATURE  AT  THE  EXACT  POSITION  OF  THE  SENSORS. 
C 

c 

DIMENSION  U(N1 ,NT)  ,X(KS) ,TEMP(KS) ,ARG( 6),VAL(6) 
DO  80  IA=1,KS 
JA=1 

AP  =  X( IA)*N+1 .0 
IAR=AR 
BR=AP-IAR 

IFCPR.GE.0.5)  GO  tq  60 
ARG( JA )  =  ( IAP-1 )*H 
VAL(JA)=U( IAR, Jl ) 
JAUX=0 
GO  TO  ~*Q 
60  ARG(JA)=TAP*H 

VAL( JA)=U( IAR+1, Jl ) 
JAUX=1 

70  CONTINUE 

IF  (JAUX.EQ.O)  GO  TO  73 
D0  71JA=2,6,2 
ARG( JA)  =  ( IAB-JA/2)*H 
VAL( JA)=U( IAR-JA/2+1, Jl) 

71  CONTINUE 

DO  12    JA=3,5,? 

ARG( JA)=(IAR+( JA-1 )/2)*H 

VAL(JA)=U(IAP+(JA-l)/2+l,Jl) 

72  CONTINUE 
GO  TO  7R 

73  DO  74  JA=2,6,2 

ARG( JA)=(  IAR  +  JA/2-la)*H 
VAL( JA)=U( IAR+JA/2, Jl) 

74  CONTINUE 

DO  ^5  JA=3,5,2 

ARG( JA)=( IAR-( JA+1)/2)*H 

VAL( JA)=U< IAP-(JA+l)/2+l,Jl) 

75  CONTINUE 

CALL  ALI  (X( IA) , AP0,VAL,TEMP( IA) ,6, .1E-03, IF°2) 
WRITE  (6,9000)  IER2 
80  CONTINUE 
9000  FORMAT  ( ' 0' , ' I ER2= ' , I  5 ) 
RETURN 
END 
C 

c 

SUBROUTINE  FEED  ( TEMP , ALPHA, KF ,KS , MP , P , Q , PHI , TAU, MI U , 
1  OMEGA, RO) 
C 

c 

C      THIS  S/R  WORKS  THE  INFORMATION  FROM  THE  SENSORS  AND 

C      COMPUTES  THE  FEEDBACK  CONTROL  LAW. 

C 

c 

REAL*4  MIU 

DIMENSION  G(MP ,KF)  ,  ALPHA  (MP)  ,  MIU(KF)  ,  P  (  KF  ,  KS  )  ,  TEM  P  (  KS  ) 
1,0MFGA(KS) ,PO(KF),°HI(KF),TAU(KF) 
CALL  MPRD  (P, TEMP, OMEGA, KF, KS, 0, 0, 1) 
DO  10  1=1, KF 
P0( I )  =  PHI(  I )*OMEGA< I  ) 
10  MIU( I)=TAU(I >-D0(I ) 

CALL  MPRD  (0, MIU, AL°HA, MP, KF  ,0,0,1) 

RETURN 

END 
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NO 


Plot  "TEMPI"  Against     "  TEMP2  " 


NO 


Compute  a 


Flow  Graph  for  the  Solution 
of  the  P.D.E.  Using  Bubnov-Galerkin  Method 
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c 

c 

c 
c 
c 


SOLUTION  BY  BU8NCV-GALERKIN  MFTHOP 


10 
11 


30 


C 
C 

c 
c 

c 

c 


REAL 
REAL 
DIME 

1P(3, 

2TEMP 
DATA 
DATA 
DATA 
1  .60, 
24.71 
3-61. 
H=2. 
N=3 
KS  =  3 
KF=3 
MIU( 
MTU( 
M I U  ( 
DO  5 
DO  5 
D=A( 
UK  I 
DO  6 
E=6. 
TFMP 
DO  1 
READ 
DO  1 
READ 
READ 
READ 
READ 
WRIT 

L  =  0 
Z10= 

Z20= 
Z30  = 
Z100 
Z200 
Z300 
CONT 
CALL 
Z100 
Z10= 
CALL 
Z200 
Z20  = 
CALL 
Z300 
Z30= 
L  =  L  + 


LAB 
*8  T 
NSIO 

3)  ,0 

1(21 

TAU 

PHI 

X/0 

.6^, 

239, 

6?50 

5   E 


1)=T 

2)=T 

3)=T 

J  =  l 

1=1 

J)*X 

,J)  = 

J  =  l 

28?1 

2(J) 

0  1  = 
(5, 

1  1  = 
(5,2 
(5,  4 
(5,4 
(5,  A 
E  (6 


ELI/4 
TITL1 
N  TAU 
MEGA( 
)  ,TE^ 
/  1 !  6  . 
/IIP. 
.0,  .0 
.70,  . 
7.8  53 
3/ 
-05 


HWANT/,LA3EL2/4HHAVE/,MIU 
(12)f  ITITL2(12)f  ITITL3U2) 

( 3)  ,Z(3) , PHI( 3)  ,TEMP (3)  ,U(3,3)  , 

3  )  , R0(3),MIU(  3)  ,P(3,3  ),  UK  2  1,3)  , 

P2(21),A(3),X(21),Y(21),A1(3) 

l,82.^,-°3.7/ 

4353, 98.6961, 59. 2177/ 

5,. 1,. 15,. 20,. 2 5,. 30,. 35,. 40,.  45,. 50,. 55, 

7  5, .80, .85, .9  0, .95, 1.0/, A/ 1 . 5  7080, 

98/,MP/21/,Al/-2.4674-0,-2  2.2  06  61, 


) 


AU(1 

AU(2) 

AU  (  3  ) 

,N 

,NP 

(T  ) 

1.41A21*SIN(D) 

,N^ 

35-X 

=  1.0- 

1,KF 

2000)     (P( I ,J) , J=1,KS) 

l.KF 

000) 

000) 

000) 

000) 

,500' 


(J ) 

>-COS(E) 


(U(  I, J )  ,J  =  1,KS) 
TTITL1 
t  TT  T|_  ^  - 
ITITL3 
'0) 


P. 

0. 

0. 

=0. 

=  0. 

=  0. 

INUE 

DIEEQ(  Z10,Z10  0,Z(  1)  ,MIU(1)  ,AK1)  ,L,H) 
=  Z10 
Z(l  ) 

DIFEQ(Z2  0,Z20  0,Z(2),MIU(2),A1(2),L,H) 
=  Z20 
Z(2) 

DIFEQ( Z30,Z300,Z( 3),MIU(3),A1(3),L,H) 
=  Z30 
Z(3) 
1 


IF  P  AMD  U  APE  SQUARE  MATRICES   P*U=I  &N0  INSTEAD  OF 
OMEGA  WE  CAN  JSE  DIRECTLY  Z. HOWEVER  A  SIMULATION  WHERE 
ERRORS  IN  THE  MEASUREMENTS  ARE  TAKEN  IN  ACCOUNT  MUST 
INCLUDE  ALL  THE  FOLLOWING  STEPS. 

50  CALL  MPRD  ( U , Z , T EMP , KS , N , 0 , 0 , 1 ) 

CALL  MPPD(U1,Z,TEM31,21,3,0,0,1) 

CALL  MPRD  (d, TEMP, OMEGA, KF,KS, 0,0, 1) 

R=L/5.0 

R 1  =  I  F I  X  (  R  ) 

R2  =  P  — R  1 

IF(R2 .LT. . 1E-04)GO  ^o  55 

GO  TO  56 
55  WRITE  (6,6000)  L , TEMPI ( 1 ) , TEMP l ( 5 ) ,TEMP! ( 9  ) , TEMPI ( 13 ) , 
1TEMPK  17)  ,TEMP1(21  ) 
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56  DO  60  1=1, KF 

RO(I  )  =  PHI ( I )*0"EGA(  !  ) 
MIU(  I)=TAU(  !)-DG(I  ) 
60  CONTINUE 

IF  ( L.EQ.100)  GO  TO  70 
IF  (L.FQ.200)  GO  TO  80 
IF  (L.FQ. 4-000)  GO  TO  90 
GO  TO  ?0 
70  CALL  DPAW( 21  ,X ,TEMP2,1 ,0,LABEL1 , IT ITL 1 ,0 , 0  , 1  ,1,2,2, 
17,7,1, LAST) 

CALL  DP. AW (21  ,X,TEMP1,3,0,LABEL2,ITITL1,0,0,1,1,2,2, 
17,7,1, LAST) 
GO  TO  30 
80  CALL  DP AW(?1, X, TEMP2, 1 ,0,LABEL1,ITITL?, 0,0,1,  1,2,2, 
17,7,1, LA  ST) 

CALL  OR AW(21 , X , TEMPI , 3 ,0 , LABFL2 , ITITL2,0,0, 1, 1,2,2, 
17,7,1, LAST) 
GO  TO  3  0 
90  CALL  DRAV'(21  ,X,TEMP2, 1,0,  LABEL1,  ITITL3,  0,0,1,  1,2,2, 
17,7,1, LAST) 

CALL  DP  AW (21, X, T EM PI, 3, 0, LAB  EL  2, ITITL 3, 0,0, 1,1, 2, 2, 
17, 7,1, LAST) 
GO  TO  100 
2000  FORMAT  (8E10.5) 
4000  FORMAT (1  GAP) 

5000  FORMAT ( « 1' ,3X,  'L' ,  9X,  ' TEMP  1 ( 1 ) ■ , 14X , ■ TF WP 1 ( 5 ) • ,  8X, 
1  "TRIP1  (9)'  ,  AX,  'TEMPK13)  '  ,  pX,«TEM.Pl(17)  '  ,  8X, 
2'TEN,PK21)  '  ,//  ) 
6000  FORMAT  ('  •  ,  I4,E?3.5,5E17.5) 
100  CONTINUE 
END 

SUBPCUTINE  DIFEQ(ZO,ZOO,Z,FIU,A,L,H) 
C 
C 

C         THIS  S/P  INTEGRATES   AN  ORDINARY  1ST  GPDEP  nIFcE_ 
C      PENTIAL  EQUATION  BY  THE  PREDICTION  AND  CORRECTION 
C      METHOD 
C 

c 

ZDOT0=A*Z00+FIU 

IF  (L.EQ.O)  ZDOT0=0. 
ZDOT1=A*ZO+F IU 
1  =  0 

L1=L+1 
10  Z1=Z00+2.0*H*ZD0T1 
1  =  1  +  1 

ZD0T1=A*Z1+FIU 
Z=ZO+H/2.0*(ZDOT1+ZDOTO) 
IF(ABS(Z-Z1) .LE.l .c-05)  GO  TO  40 
IF( I .ED. 100)  GO  TO  3^ 
GO  TO  10 

39  WRITE  (6,1000)  LI 

40  RETURN 

1000    FORMAT     ('0','NO    CONVERGENCE «, 3X  ,' L 1  =  ', 1 4  ,/  ) 
END 
//GO.SYSIN    DD    * 

.17532  EOO. 31593  EOO. 393^1  EOO 
.39393  EOO. 17533  E00-. 31 5°1 E 00 
.31589  F00-.393O1E00. 1753?  EOO 
.61357  EOC. 13787  E01. 11058  E01 
.11057  E01. 61353  EOO- . 13 7R«F01 
.13788     F01-.11057E01. 61367    EOO 

WANTED    AND    OBTAINED    TEMPERATURE    AT    0.0025    SEC  ALM0403    A. MO 

-BUBNOV-GALERKIN 

WANTED    AND    OBTAINED    TEMPERATURE    AT    0.0050    SEC  ALM0403    A. MO 

-BUBNOV-GALEPKIN 

WAN'TED    AND    CRTAINEO    TEMPERATURE    AT    0.1  SEC  ALM0403    A. MO 

-BUBNOV-GALEPKIN 


159 


c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 
c 
c 
c 
c 


THIS  P 

SOLVES 

METHOD 

OF  THE 

PLOT  0 

OPTIMI 

OBTAIN 

MODEL. 

D I  MENS 

COMMON 

COMMON 

EXTERN 

NW=3 

TAU( 1) 

TAU(2) 

TAU(3) 

X(l)=. 

X(2)=. 

X(3)=. 

X(4)  =  . 

X(5)=. 

X(6)=. 


0(3,3) ,PHI (3) ,X(6) 


ROGRAM  IS  BASICALLY  THE  SAME  AS  THE  ONE  THAT 

THE  CL0SE-LOOD  PROBLEM  USING  THE  CR ANK-NI COL  SON 
.THE  MAIN  DIFFERENCE  IS  IN  THE  FACT  THAT  MOST 

WRITE  STATEMENTS  WERE  REMOVED, AS  WELL  AS  THE 
UTPUTS.BY  DEFINING  A  COST  FUNCTION  WHICH  IS 
ZED  BY  A  GRADIENT  SEARCH  A  SET  OP  TAL'S  IS 
ED  AMD  USED  AS  INPUT  TO  THE  SYSTEM  IN  THE  BASIC 
ONLY  THEM  A  PLOT  OUTPUT  IS  OBTAINED 
ION  T£U(3) 
/0NE/P(3,6) , 
/TWO/N,M,KF,KS,MD,NT 
AL  SQUERR 

=116.9 

=  38.2 

=-35.6 

1540  6 

30828 

46234 

61621 

76996 

92325 


N 

M 

TION 
KF  — 
KS  — 
MP- 
N  =  20 
M=4 
KF  =  3 
KS=6 
MP=3 
NT  =  4 


-#  OF  EQUATIONS 

-#  OF  RIGHT  HAND  VECTORS  IN  THF  SYSTEM  OF  EQUA 

S  T0  BE  SOLVED  BY  CR ANK-NICOLSON  METHOD. 

— #  OF  FILTERS 

— #    OF  SENSORS 

— #     OF  MANIPULATORS 


0 


PHK  1)  =  118.4353 
PHI (2)=°8.6961 
PHI(3)=5°.2177 

DEFINE   THE  MATRICES  P  AND  Q  OF  THE  FEEDBACK  LOOP 


DO  4  I=1,KF 

READ  (5,2000) 
DO  5  1  =  1, MP 
READ  (5,2000) 
DELCAP=10. 
CALL  DIRECT(TAU, MW 
1K0NVRG, 100,-1) 
WRITE(6,300) 
WRITE(6,400)KONVRG 
WRITE(6, 500) 


(P(I,J)  ,J=1,KS) 
(0(1  ,  J)  ,  J=1,KF) 


COS  MI Nt DEL CAP, .25 » .01 , SQUERR 


10 

300 
400 
500 
600 
700 
2000 


DO  10  I=1,NW 


WRITE(6 
WRITE(6 
FORMAT ( 
FORMAT ( 
FORMAT ( 
FORMAT ( 
FORMAT ( 
FORMAT 
STOP 
END 


600) 
700) 


TAUU  ) 

COSMIN 


0' 
0» 
0' 
0' 
0« 


,T60, 
,T60, 
tT60, 
,T50, 
*T60, 


3E10.5 ) 


' KONVRG' ) 

14,/) 

•TAU' ) 

F15.5) 

• SQUERR=« 


F7.3) 


C 
C 
C 
C 
C 
C 


FUNCTION  SQUERR  (TAU) 


MODAL  CONTROL  OF  A  DISTRIBUTED  PARAMETER  SYSTEM 
DESCRIBED   BY  A  PARABOLIC  DIFFERENTIAL  EQUATION 
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c 

C         CARDS  UP  TO  "HO  4 — •'  AND  THE  DATA  CARDS  MUST  RE 

C  CHANGED  ACCORDING  TO  WHAT  IS  REQUIRED  FROM  THE 

C  SYSTEM  AND  ALSO  ACCORDING  TO  THE  PHYSICAL  LIMITA 

C         TIONS.THE  ONLY  ADDITIONAL  CHANGES  MAY  RE  IN  LOOP' 

C  "DO  38 — ", WHICH  DEFINES  THE  HEAT  DISTRIBUTION  OF 

C         THE  MANIDULATORS, IN  STATEMENT  «4?,IN  THE  LOOP 

C  "DO  41 — "  AND  IN  THE  PARAMETERS  AND  VECTOR  Y2  OF 

C         S/R  "DRAW". 

C 

C 

REAL*8  R,A1 

REAL*4  MIU 

DIMENSION  ALPHA(3) ,U(21,40   ),R(20)fA(59  ), 
1EVK21)  ,EV2(21)  ,EV3(21)  ,H1(21)  ,H2(21)f 
2H3(21) ,TAU(3),XEMP(6),X1(21),Y1(21),A1(59), 
3RK80)  , R2(20),MIU(3)  ,  OMEGA  (  3  )  ,  RO  (  3  )  ,Y2(21),SQU(21) 

,6) ,0(3,3) , 
,KF, KS,MP,NT 

N1=N+1 

N2=N+2 

H=1.0/N 

HK=H**2 

Jl=l 
C 

C  BOUNDAPY    CONDITIONS    AT    X=0 

C 

DO    2    J=1,NT 
2    U(1,J)=0.0 
C 
C  INITIAL    CONDITIONS 


C 


C 


DO    3    I=1,N1 
U(  1, 1)=0.0 


C 

C  SOLUTION    BY    CRANK-NICHOLSON    METHOD 

C 

c 

C      DEFINING  THE  TRIDIAGONAL  MATRIX  OF  EQUATIONS  ACCOR- 

C      DING  TO  THE  REQUIREMENTS  OF  S/R  DGELB. 

C 

NM=N*(M-1) 

NM0=NM-2 

NO=NM-? 

DO  30  I=1,NM0,3 

A(  I  )=4.0 

IF(I.EQ.NMO)  GO  TO  30 

11=1+1 

12=1+2 

A( Il)=-1.0 

A( I2)=-1.0 
30  CONTINUE 

A(N0)=-2.0 
C 

C      DEFINING  THE  FIRST  (M-1)*N  TERMS  OF  Rl 
C 

DO  35  IA=1,NM 
35  R1(IA)=0.0 
C 

C      HEAT  DISTRIBUTION  OF  THE  MAN  I PUL ATORS .THF  NEXT  CARDS, 
C      UP  TO  STATEMENT  #38  ^UST  BE  CHANGED  IF  THE  niSTDIPU- 
C      TION  OF  HEAT  CORRESPONDENT  TO  EACH  MANIPULATOR  ALSO 
C      CHANGES. 


DO  38  ID=1,N1 
EVK  ID)  =  1.570P0*H*(  ID-1) 
EV2UD)=4.7123  9*H*(  ID-1) 
EV3( ID) =7. 9  5398*H*( TD-1) 
H1(ID)=  1.414214*SIN(EV1 ( ID) ) 
H2(ID)=  1.414214*SIN(EV?(  ID)  ) 
38  H3(ID)=  1.414214*SIN(EV3( ID) ) 
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c 
c 
c 


c 
c 
c 
c 


c 
c 

c 


c 
c 
c 
c 


IF  I.C.'S.NE.O  INSTEAD  OF  TAU  WE  MUST  COMPUTE  AND  USE 
MIU. 

CALL  MPRD(Q,TAU, AL PHA , MP , KF , 0 , 0 , 1 ) 
40  CONTINUE 

42  DO  43  1  =  1, N 
ID=I+1 
IE=NM+I 

43  R1(IF)=( ALPHA! 1) *H1 ( I D ) + ALPHA ( 2 ) *H2 ( ID) +ALPHA( 3) * 
1H3UD)  )*HK*2.0 

S/R  DGELB  SOLVES  IN  DOUBLE  PRECISION  A  SYSTEM  OF 
LINEAR  EQUATIONS. 

DH  44  1=1, NMO 

44  Aid  )=A(  I) 
N3=2*N+1 
N4=3*N+1 

DO  45  I=1,M 

R2(I  )  =  R1(  I)+R1 (M1)+R1(N3)+R1(N4) 

N1=N1+1 

N3=N3+1 

45  N4=N4+1 

DO  47  I=1,N 
R( I)=P2(  I) 

47  CONTINUE 
MUD=1 
MLD=1 

CALL    DGELB    ( P , A  1 , N, 1 , MUD ,MLD, . 1 E-06, I ER 1 ) 

REDEFINING    THE     I.C.'S 

DO    48    ID=2,N 

48  PI  (ID) =0.0 
NE=2*N-1 
N0=N-1 
NE1=NE+1 
IG=1 

DO    49     ID=N2,NE 
R1(ID)=R( IG) 

49  IG=IG+1 
R1(NE1)=2*R(N0) 
NF=3*N 
NE2=NE+2 

IG=2 

DO    50    ID=NP2,NF 

R1(ID)=R( IG) 

50  IG=IG+1 
R1(NF)=0.0 
N1=N+1 

DO  60  ID=2,N1 
IE=ID-1 
60  U(ID,J1)=R(  IE) 

COMPUTATION  OF  VECTORS  ARG  AND  VAL  FOR  THE  INTERPC. 
LATING  S/R  "ALI". 

CALL  ANA  (TEMP,KS,J1,H,N,U,X,N1,NT) 

CALL  FEED  (  TEMP, ALPHA, KF, KS, MP, P,Q, PHI , TAU, MIU, OMEGA, 
1R0) 
J1=J1+1 

IFU1.EQ.38)  WRITE(6,6000)  J  1  ,  (  ALPHA(  I  )  ,  1=  1 ,  MP  ) 
IF(Jl.GT.NT)  30  TO  90 
GO  TO  40 
90  IQ=40 

WRITE  (6,6000)  J  1 , ( AL PHA ( I ) , I  =  1 , MP) 

SQUERR=0. 

DO  9  5  IH=1,N1 

X1(IH)=(IH-1)*H 

D=6.283185-X1(  IH) 

Y2( IH)=1.0-C0S (D) 

Yl(  IH)=U(IH, IQ) 
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SQU( IH)  =  (Y1(  IH)-Y2( IH)  )**2 
95  SPUERR=SQUERR+SQU( IH) 
6000  FORMAT  ( ■ 0 ■ , ' ALPHA ( • , 16 , • ) = • , 7E 14. 5) 
.  END 
C 
C 

SUBROUTINE  ANA  ( TEMP, KS , Jl , H , N, U , X »N1 ,NT ) 
C 

c 

C  THIS  S/P  COMPUTES  THE  INTERPOLATED  VALUES  OF  THE 
C  TEMPERATURE  AT  THE  EXACT  POSITION  OF  THE  SENSORS 
C 

c 

DIMENSION  U(N1,NT),X(KS)  ,TEMP(KS)  ,ARG(6) ,VAL(6) 
DO  80  IA=1,KS 
JA  =  1 

AR=X(I£ )*N+1.0 
IAR=AR 
BR=AR-IAR 

IF(BR.OE.0.5)  Gn  TO  60 
ARG( JA )=(IAR-1 )*H 
VAL( JA)=U( IAP, Jl ) 
JAUX^O 
GO  TO  70 
60  ARG( JA)=IAP-H 

VAL(JA)=U(  IAR+1, Jl  ) 
JAUX=1 

70  CONTINUE 

IF  (JAUX.EO.O)  GO  TO  73 

DO  71  JA=2,6,2 

ARG( JA)=( IAR-JA/2) *H 

VAL(JA)=U(IAR-JA/2+l,Jl) 

71  CONTINUE 

DH  72  JA=3,5,2 
ARG(JA)=(IAP+( JA-1 )/?)*H 
VAL( JA)=U(  IAR-K  JA-1  )/2  +  l,Jl) 

72  CONTINUE 
GO  TO  75 

73  DO  74  JA=2,6,2 
ARG(JA)=(I AR+JA/2-l)*H 
VAL(  JA)=U(  IAP+JA/2, Jl  ) 

74  CONTINUE 

DO  75  JA=3,5,2 

APGC JA)=(IAR-( JA+1 )/2)*H 

VAL(  JA)=U(  IAR-(  JA+D/2+1,  Jl) 

75  CONTINUE 

CALL  ALI  (X(IA)  , ARG, VAL , TEMP ( I  A) ,6, .1E-03,  IER2) 
PO  CONTINUE 

RETURN 

END 
C 
C 


SUBRHUTINE  FEED  (TEMP t ALPHA, KF , KS , MP t P , 0 , PHI , TAU , MIU, 
IOMEGA, RO) 
C 
C 

C      THIS  S/R  WORKS  THE  INFORMATION  fp.CM  THE  SENSORS  AND 
C      COMPUTES  THE  FEEDBACK  CONTROL  LAW. 
C 

c 

REAL*4  MIU 

DIMENSION  Q(MPtKF) , ALPHA (MP) ,MIU(KF),P(KF,KS),TEMD(KS) 
1,0MEGA( KS) ,P0( KF) ,°HI (Kc) ,TAU(KF) 
CALL  MPPD  (P, TEM°, OMEGA, KF, KS, 0, 0, 1 ) 
DO  10  1=1, KF 
R0( I )=PHI( I )*OMFGA( I ) 
10  MIU(  I)=TAU(  I )-PO(I  ) 

CALL  MPRD  (0, MIU, ALPHA, MP, KF, 0,0,1) 
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